#############################################################
# Naima Green-Riley and Andrew Leber
# Whose War is it Anyway?
# Article about racial differences in support for war
# Last edited: November 23, 2022
#
# 1) Load required data
# 2) Race-gap regressions & figures
# 3) Over-time plots 
# 4) Race gap (Democrats only)
# 5) Risks of War (Casualties & Military Service)
# 6) Costs of War (Spending)
# 7) War & Alienation
# 8) Plot of all explanatory variables
# 9) Afghanistan replication, overall summary plots
# 10) Afghanistan replication, explanatory variables
# 11) Robustness check: No income control variable
# 12) Robustness check: Alternative Casualties Measure
# 13) Robustness check: Risks of War, bivariate
# 14) Robustness check: Costs of War, bivariate
# 15) Robustness check: Alienation and War, bivariate
# 16) Mediations sensitivity analysis
# 17) Specific statistics (various)
# 18) Gulf War barplot
# 
#############################################################

## Quick fix for stargazer <= 5.2.3 is.na() issue with long model names in R >= 4.2
# Solution by Alexey Knorre: https://gist.github.com/alexeyknorre/b0780836f4cec04d41a863a683f91b53

# Unload stargazer if loaded
detach("package:stargazer",unload=T)
# Delete it
remove.packages("stargazer")
# Download the source
download.file("https://cran.r-project.org/src/contrib/stargazer_5.2.3.tar.gz", destfile = "stargazer_5.2.3.tar.gz")
# Unpack
untar("stargazer_5.2.3.tar.gz")
# Read the sourcefile with .inside.bracket fun
stargazer_src <- readLines("stargazer/R/stargazer-internal.R")
# Move the length check 5 lines up so it precedes is.na(.)
stargazer_src[1990] <- stargazer_src[1995]
stargazer_src[1995] <- ""
# Save back
writeLines(stargazer_src, con="stargazer/R/stargazer-internal.R")
# Compile and install the patched package
install.packages("stargazer", repos = NULL, type="source")

# Required packages
attach(mtcars)
require(foreign)
library(memisc)
library(mediation)
library(foreign)
library(ggplot2)
library(car)
library(stargazer)
library(lmtest)
library(vtable)
library(showtext)

#############################################################
# 1) Read in datasets for article
#############################################################

# Set working directory to wherever downloaded
# ENTER WORKING DIRECTORY

# CCES data (military spending & military families)
cces_bush<-read.csv("cces_bush2_new.csv", check.names = F, stringsAsFactors = F)
cces_obama<-read.csv("cces_obama2_new.csv", check.names = F, stringsAsFactors = F)

# Asstd. polls (alienation, Afghanistan)
Bush_data_BW<-read.csv("alienation_bush_afghan_new.csv", check.names = F, stringsAsFactors = F)
Bush_data_BW<-Bush_data_BW[which(Bush_data_BW$black%in%c(0,1)),]

Obama_data_BW<-read.csv("alienation_obama_afghan_new.csv", check.names = F, stringsAsFactors = F)
Obama_data_BW<-Obama_data_BW[which(Obama_data_BW$black%in%c(0,1)),]

# Asstd. polls (alienation, Iraq)
Bush_data_IQ_BW<-read.csv("alienation_bush_iraq_new.csv", check.names = F, stringsAsFactors = F)
Bush_data_IQ_BW<-Bush_data_IQ_BW[which(Bush_data_IQ_BW$black%in%c(0,1)),]
Obama_data_IQ_BW<-read.csv("alienation_obama_iraq_new.csv", check.names = F, stringsAsFactors = F)
Obama_data_IQ_BW<-Obama_data_IQ_BW[which(Obama_data_IQ_BW$black%in%c(0,1)),]

# Oct 2005 CNN poll (Bush, Iraq)
Bush_IQ_CNN2005 <- read.csv("Bush_IQ_CNN2005.csv", check.names = F, stringsAsFactors = F)
  
# ABC polls (casualties)
abc_iraq_combine<-read.csv("abc_casualties_iraq2_new.csv", check.names = F, stringsAsFactors = F)
abc_iraq_combine<-abc_iraq_combine[which(abc_iraq_combine$black%in%c(0,1)),]
poll_afgh<-read.csv("abc_casualties_afgh2_new.csv", check.names = F, stringsAsFactors = F)
poll_afgh<-poll_afgh[which(poll_afgh$black%in%c(0,1)),]

# standardize poll FE variable
abc_iraq_combine$poll<-factor(abc_iraq_combine$poll_code)
cces_bush$poll<-factor(cces_bush$year)
cces_obama$poll<-factor(cces_obama$year)
Bush_data_BW$poll<-factor(Bush_data_BW$poll)
Bush_data_IQ_BW$poll<-factor(Bush_data_IQ_BW$poll)
Obama_data_BW$poll<-factor(Obama_data_BW$poll)
Obama_data_IQ_BW$poll<-factor(Obama_data_IQ_BW$poll)


# Relevel factors
cces_bush$income<-relevel(factor(cces_bush$income), ref = "Less than $30K")
cces_obama$income<-relevel(factor(cces_obama$income), ref = "Less than $30K")
Bush_data_BW$income<-relevel(factor(Bush_data_BW$income), ref = "Less than $30K")
Bush_data_IQ_BW$income<-relevel(factor(Bush_data_IQ_BW$income), ref = "Less than $30K")
Obama_data_BW$income<-relevel(factor(Obama_data_BW$income), ref = "Less than $30K")
Obama_data_IQ_BW$income<-relevel(factor(Obama_data_IQ_BW$income), ref = "Less than $30K")
abc_iraq_combine$income<-relevel(factor(abc_iraq_combine$income), ref = "Less than $30K")
Bush_IQ_CNN2005$income<-relevel(factor(Bush_IQ_CNN2005$income), ref = "Less than $30K")

cces_bush$education<-factor(cces_bush$education,
                            levels = c("HS or Less", "Some College", "College Grad", "Post-Grad"))
cces_obama$education<-factor(cces_obama$education,
                                     levels = c("HS or Less", "Some College", "College Grad", "Post-Grad"))
Bush_data_BW$education<-factor(Bush_data_BW$education,
                                       levels = c("HS or Less", "Some College", "College Grad", "Post-Grad"))
Bush_data_IQ_BW$education<-factor(Bush_data_IQ_BW$education,
                                          levels = c("HS or Less", "Some College", "College Grad", "Post-Grad"))
Obama_data_BW$education<-factor(Obama_data_BW$education,
                                        levels = c("HS or Less", "Some College", "College Grad", "Post-Grad"))
Obama_data_IQ_BW$education<-factor(Obama_data_IQ_BW$education,
                                           levels = c("HS or Less", "Some College", "College Grad", "Post-Grad"))
abc_iraq_combine$education<-factor(abc_iraq_combine$education,
                                           levels = c("HS or Less", "Some College", "College Grad", "Post-Grad"))
Bush_IQ_CNN2005$education<-factor(Bush_IQ_CNN2005$education,
                                          levels = c("HS or Less", "Some College", "College Grad", "Post-Grad"))
poll_afgh$education<-relevel(factor(poll_afgh$education), ref = "No College")


### set working directory to output for graphs

# ENTER FOLDER

##

# load font for plot output
font_add_google("EB Garamond")
showtext_auto()
 

####################################################
# 1a) Generate Summary tables for Datasets
####################################################

# Casualty
sumtable(abc_iraq_combine, vars = c("support_war", "casual_sens", "casualty_unaccept",
                                    "black","gender", 
                                      "education","party_id","age","income"),
         out='latex',file='1_abc_iraq_casualty_sum.tex',factor.numeric = TRUE, digits = 2, fixed.digits = TRUE,
         summ=c('mean(x)', 'sd(x)','min(x)', 'max(x)'))

sumtable(poll_afgh, vars = c("support_war", "casualty_sens", "black","gender", 
                                    "education","party_id","age"),
         out='latex',file='2_abc_afg_casualty_sum.tex',factor.numeric = TRUE, digits = 2, fixed.digits = TRUE,
         summ=c('mean(x)', 'sd(x)','min(x)', 'max(x)'))

# Iraq (Alienation) - Asstd.
sumtable(Bush_data_IQ_BW, vars = c("outcome", "alienation_binary", "black","gender", 
                                    "education","party_id","age","income"),
         out='latex',file='3_asstd_iraq_bush_alien_sum.tex',factor.numeric = TRUE, digits = 2, fixed.digits = TRUE,
         summ=c('mean(x)', 'sd(x)','min(x)', 'max(x)'))

sumtable(Obama_data_IQ_BW, vars = c("outcome", "alienation_binary", "black","gender", 
                                    "education","party_id","age","income"),
         out='latex',file='4_asstd_iraq_obama_alienation_sum.tex',factor.numeric = TRUE, digits = 2, fixed.digits = TRUE,
         summ=c('mean(x)', 'sd(x)','min(x)', 'max(x)'))

# Afghanistan (Alienation) - Asstd.
sumtable(Bush_data_BW, vars = c("outcome", "alienation_binary", "black","gender", 
                                   "education","party_id","age","income"),
         out='latex',file='5_asstd_afgh_bush_alien_sum.tex',factor.numeric = TRUE, digits = 2, fixed.digits = TRUE,
         summ=c('mean(x)', 'sd(x)','min(x)', 'max(x)'))

sumtable(Obama_data_BW, vars = c("outcome", "alienation_binary", "black","gender", 
                                    "education","party_id","age","income"),
         out='latex',file='6_asstd_afgh_obama_alien_sum.tex',factor.numeric = TRUE, digits = 2, fixed.digits = TRUE,
         summ=c('mean(x)', 'sd(x)','min(x)', 'max(x)'))

# CCES - Iraq - Bush
sumtable(cces_bush, vars = c("iraq_support", "serving", "dom_over_def",
                             "black","gender", 
                                "education","party_id","age","income"),
         out='latex',file='7_cces_iraq_bush_sum.tex',factor.numeric = TRUE, digits = 2, fixed.digits = TRUE,
         summ=c('mean(x)', 'sd(x)','min(x)', 'max(x)'))

# CCES - Iraq - Obama
sumtable(cces_obama, vars = c("iraq_support", "serving", "dom_over_def",
                             "black","gender", 
                             "education","party_id","age","income"),
         out='latex',file='8_cces_iraq_obama_sum.tex',factor.numeric = TRUE, digits = 2, fixed.digits = TRUE,
         summ=c('mean(x)', 'sd(x)','min(x)', 'max(x)'))

# CCES - Afghanistan - Obama
sumtable(cces_obama, 
         vars = c("afghanistan_support", "serving", "dom_over_def",
                              "black","gender", 
                              "education","party_id","age","income"),
         out='latex',file='9_cces_afgh_obama_sum.tex',factor.numeric = TRUE, digits = 2, fixed.digits = TRUE,
         summ=c('mean(x)', 'sd(x)','min(x)', 'max(x)'))

###############################################################
# 2) Regressions for Table 2, Figure 3 (Race Gap, Iraq War)
###############################################################

# Models 1 & 2: Between-group & with controls, ABC polling dataset
bush_iraq_casualty_reg_1<-lm(support_war ~ black+factor(poll), 
                             data = abc_iraq_combine)

bush_iraq_casualty_se_1<-sqrt(diag(vcovHC(bush_iraq_casualty_reg_1, type = "HC1")))

bush_iraq_casualty_reg_2<-lm(support_war ~ black+gender+party_id+education+age+income+
                                         factor(poll), 
                                       data = abc_iraq_combine)

bush_iraq_casualty_se_2<-sqrt(diag(vcovHC(bush_iraq_casualty_reg_2, type = "HC1")))

# Models 3 & 4: Between-group & with controls, Asstd. (Alienation) polling dataset, Bush
bush_iraq_reg_3<-lm(outcome ~ black+factor(poll), 
                    data = Bush_data_IQ_BW)
bush_iraq_se_3<-sqrt(diag(vcovHC(bush_iraq_reg_3, type = "HC1")))

bush_iraq_reg_4<-lm(outcome ~ black+gender+party_id+education+age+income+factor(poll), 
                    data = Bush_data_IQ_BW)
bush_iraq_se_4<-sqrt(diag(vcovHC(bush_iraq_reg_4, type = "HC1")))

# Models 5 & 6: Between-group & with controls, Asstd. (Alienation) polling dataset, Obama
obama_iraq_reg_5<-lm(outcome ~ black+factor(poll), 
                     data = Obama_data_IQ_BW)
obama_iraq_se_5<-sqrt(diag(vcovHC(obama_iraq_reg_5, type = "HC1")))

obama_iraq_reg_6<-lm(outcome ~ black+gender+party_id+education+age+income+factor(poll), 
                     data = Obama_data_IQ_BW)
obama_iraq_se_6<-sqrt(diag(vcovHC(obama_iraq_reg_6, type = "HC1")))

# Models 7 & 8: Between-group & with controls, CCES polling dataset, Bush
iraq_bush_cces_reg_7<-lm(iraq_support ~ black + factor(poll), data = cces_bush)

iraq_bush_cces_se_7<-sqrt(diag(vcovHC(iraq_bush_cces_reg_7,type = "HC1")))

iraq_bush_cces_reg_8<-lm(iraq_support ~ black+gender+party_id+education+age+income+
                           factor(poll), data = cces_bush)

iraq_bush_cces_se_8<-sqrt(diag(vcovHC(iraq_bush_cces_reg_8,type = "HC1")))

# Models 9 & 10: Between-group & with controls, CCES polling dataset, Obama
iraq_obama_cces_reg_9<-lm(iraq_support ~ black + factor(poll), data = cces_obama)

iraq_obama_cces_se_9<-sqrt(diag(vcovHC(iraq_obama_cces_reg_9,type = "HC1")))

iraq_obama_cces_reg_10<-lm(iraq_support ~ black+gender+party_id+education+age+income+
                             factor(poll), data = cces_obama)

iraq_obama_cces_se_10<-sqrt(diag(vcovHC(iraq_obama_cces_reg_10,type = "HC1")))

####################################################
# 2a) Output for Table 2 (Regression Models 1-10)
####################################################
stargazer(bush_iraq_casualty_reg_1, bush_iraq_casualty_reg_2, 
          bush_iraq_reg_3, bush_iraq_reg_4, 
          obama_iraq_reg_5, obama_iraq_reg_6, 
          iraq_bush_cces_reg_7, iraq_bush_cces_reg_8,
          iraq_obama_cces_reg_9, iraq_obama_cces_reg_10,
          se = list(bush_iraq_casualty_se_1, bush_iraq_casualty_se_2, 
                    bush_iraq_se_3, bush_iraq_se_4, 
                    obama_iraq_se_5, obama_iraq_se_6,
                    iraq_bush_cces_se_7, iraq_bush_cces_se_8, 
                    iraq_obama_cces_se_9, iraq_obama_cces_se_10),
          covariate.labels = c("Black", "Male",
                               "Independent", "Republican",
                               "Some College", "College Grad", "Post-Grad",
                               "Age: 30-49", "Age: 50-64", "Age: 65+",
                               "Income: $30K-$50K",
                               "Income: $50K-$75K",
                               "Income: $75K and above",
                               rep(NA, 61), "Constant"),
          omit = c("poll"),
          omit.labels = c("Poll FE"),
          column.labels = c("Bush (Casualties)", "Bush (Alienation)",
                            "Obama (Alienation)", "Bush (CCES)" ,"Obama (CCES)"),
          column.separate = c(2,2,2,2,2),
          add.lines = list(c("Polls", 18, 17, 39, 35,
                             11,11,2,2,3,3)),
          omit.stat = c("adj.rsq","f", "ser"),
          out = "table_2_full_regression_new.tex")


####################################################
# 2b) Output for Figure 3 (Coefficient Plot, "Black," Models 1-10)
####################################################

# Coefficients from each model for "black"
coef_black<-c(coef(bush_iraq_casualty_reg_1)[2],
              coef(bush_iraq_casualty_reg_2)[2],
              coef(bush_iraq_reg_3)[2],
              coef(bush_iraq_reg_4)[2],  
              coef(obama_iraq_reg_5)[2], 
              coef(obama_iraq_reg_6)[2], 
              coef(iraq_bush_cces_reg_7)[2],
              coef(iraq_bush_cces_reg_8)[2], 
              coef(iraq_obama_cces_reg_9)[2],
              coef(iraq_obama_cces_reg_10)[2])

# SE from each model for "black"
se_black<-c(bush_iraq_casualty_se_1[2],
            bush_iraq_casualty_se_2[2], 
            bush_iraq_se_3[2],
            bush_iraq_se_4[2], 
            obama_iraq_se_5[2],
            obama_iraq_se_6[2], 
            iraq_bush_cces_se_7[2],
            iraq_bush_cces_se_8[2], 
            iraq_obama_cces_se_9[2],
            iraq_obama_cces_se_10[2])

# Labels for each dataset used
poll_names<-c("ABC\n(Bush, Iraq)", "ABC\n(Bush, Iraq)",
              "Asstd.\n(Bush, Iraq)", "Asstd.\n(Bush, Iraq)", 
              "Asstd.\n(Obama, Iraq)", "Asstd.\n(Obama, Iraq)",
              "CCES\n(Bush, Iraq)", "CCES\n(Bush, Iraq)",
              "CCES\n(Obama, Iraq)","CCES\n(Obama, Iraq)")

# Identifiers
poll_type<-c("Bivariate", "Full Regression",
             "Bivariate", "Full Regression",
             "Bivariate", "Full Regression",
             "Bivariate", "Full Regression",
             "Bivariate", "Full Regression")

# Combine into data frame for visualization
dat<-data.frame(names = poll_names,
                type = poll_type,
             coef = coef_black,
             se = se_black)

# Generate upper and lower bounds for error bars
dat$ub<-dat$coef+1.96*dat$se
dat$lb<-dat$coef-1.96*dat$se

dat$ub90<-dat$coef+1.645*dat$se
dat$lb90<-dat$coef-1.645*dat$se

# Order the polls
dat$names=factor(dat$names, levels = c("ABC\n(Bush, Iraq)", 
                                       "Asstd.\n(Bush, Iraq)", 
                                       "Asstd.\n(Obama, Iraq)", 
                                       "CCES\n(Bush, Iraq)", "CCES\n(Obama, Iraq)"))

# Generate graph for publication
pdf(file="black_coef_new_garamond.pdf", width = 10, height = 5)
ggplot(dat, aes(names, coef, group = type, color = type))+
  geom_errorbar(ymin = dat$lb, ymax = dat$ub, width = 0.2, 
                position=position_dodge(width=0.5))+
  geom_point(position=position_dodge(width=0.5))+
  scale_y_continuous(name = "Coefficient on 'Black'",limits = c(-0.4, 0.2),
                     breaks = c(-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2))+
  geom_hline(yintercept = 0)+scale_x_discrete(name = "")+
  scale_color_manual(name = "Specification", values = c("grey55", "black"))+theme_bw()+
  theme_set(theme_bw(base_size = 20, base_family = "EB Garamond" )) 
dev.off()

setEPS()                                             
postscript("black_coef_new_garamond.eps", width = 10, height = 5) 
ggplot(dat, aes(names, coef, group = type, color = type))+
  geom_errorbar(ymin = dat$lb, ymax = dat$ub, width = 0.2, 
                position=position_dodge(width=0.5))+
  geom_point(position=position_dodge(width=0.5))+
  scale_y_continuous(name = "Coefficient on 'Black'",limits = c(-0.4, 0.2),
                     breaks = c(-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2))+
  geom_hline(yintercept = 0)+scale_x_discrete(name = "")+
  scale_color_manual(name = "Specification", values = c("grey55", "black"))+theme_bw()+
  theme_set(theme_bw(base_size = 20, base_family = "EB Garamond" )) 
dev.off()


######################################################################
# 2c) Code for Table 22, Figure 11 - As above, Democrats only
######################################################################

# Models 1 & 2: Between-group & with controls, ABC polling dataset
bush_iraq_casualty_dem_reg_1<-lm(support_war ~ black+factor(poll), 
                              data = abc_iraq_combine[which(abc_iraq_combine$party_id=="democrat"),])

bush_iraq_casualty_dem_se_1<-sqrt(diag(vcovHC(bush_iraq_casualty_dem_reg_1, type = "HC1")))

bush_iraq_casualty_dem_reg_2<-lm(support_war ~ black+gender+education+age+income+
                               factor(poll), 
                             data = abc_iraq_combine[which(abc_iraq_combine$party_id=="democrat"),])

bush_iraq_casualty_dem_se_2<-sqrt(diag(vcovHC(bush_iraq_casualty_dem_reg_2, type = "HC1")))

# Models 3 & 4: Between-group & with controls, Asstd. polling dataset, Bush
bush_iraq_dem_reg_3<-lm(outcome ~ black+factor(poll), 
                     data = Bush_data_IQ_BW[which(Bush_data_IQ_BW$party_id=="democrat"),])
bush_iraq_dem_se_3<-sqrt(diag(vcovHC(bush_iraq_dem_reg_3, type = "HC1")))

bush_iraq_dem_reg_4<-lm(outcome ~ black+gender+education+age+income+factor(poll), 
                    data = Bush_data_IQ_BW[which(Bush_data_IQ_BW$party_id=="democrat"),])
bush_iraq_dem_se_4<-sqrt(diag(vcovHC(bush_iraq_dem_reg_4, type = "HC1")))

# Models 5 & 6: Between-group & with controls, Asstd. polling dataset, Obama
obama_iraq_dem_reg_5<-lm(outcome ~ black+factor(poll), 
                      data = Obama_data_IQ_BW[which(Obama_data_IQ_BW$party_id=="democrat"),])
obama_iraq_dem_se_5<-sqrt(diag(vcovHC(obama_iraq_dem_reg_5, type = "HC1")))

obama_iraq_dem_reg_6<-lm(outcome ~ black+gender+education+age+income+factor(poll), 
                     data = Obama_data_IQ_BW[which(Obama_data_IQ_BW$party_id=="democrat"),])
obama_iraq_dem_se_6<-sqrt(diag(vcovHC(obama_iraq_dem_reg_6, type = "HC1")))

# Models 7 & 8: Between-group & with controls, CCES (Bush)
iraq_bush_cces_dem_reg_7<-lm(iraq_support ~ black + factor(poll), 
                              data = cces_bush[which(cces_bush$party_id=="democrat"),])

iraq_bush_cces_dem_se_7<-sqrt(diag(vcovHC(iraq_bush_cces_dem_reg_7,type = "HC1")))

iraq_bush_cces_dem_reg_8<-lm(iraq_support ~ black+gender+education+age+income+
                               factor(poll), data = cces_bush[which(cces_bush$party_id=="democrat"),])

iraq_bush_cces_dem_se_8<-sqrt(diag(vcovHC(iraq_bush_cces_dem_reg_8,type = "HC1")))

# Models 9 & 10: Between-group & with controls, CCES (Obama)
iraq_obama_cces_dem_reg_9<-lm(iraq_support ~ black + factor(poll), 
                               data = cces_obama[which(cces_obama$party_id=="democrat"),])

iraq_obama_cces_dem_se_9<-sqrt(diag(vcovHC(iraq_obama_cces_dem_reg_9,type = "HC1")))

iraq_obama_cces_dem_reg_10<-lm(iraq_support ~ black+gender+education+age+income+
                                 factor(poll), data = cces_obama[which(cces_obama$party_id=="democrat"),])

iraq_obama_cces_dem_se_10<-sqrt(diag(vcovHC(iraq_obama_cces_dem_reg_10,type = "HC1")))

#######################################################################
# 2d) Output for Table 22 in the Appendix (Regression Models 1-10)
#######################################################################
stargazer(bush_iraq_casualty_dem_reg_1, bush_iraq_casualty_dem_reg_2, 
          bush_iraq_dem_reg_3, bush_iraq_dem_reg_4, 
          obama_iraq_dem_reg_5, obama_iraq_dem_reg_6, 
          iraq_bush_cces_dem_reg_7, iraq_bush_cces_dem_reg_8,
          iraq_obama_cces_dem_reg_9, iraq_obama_cces_dem_reg_10,
          se = list(bush_iraq_casualty_dem_se_1, bush_iraq_casualty_dem_se_2, 
                    bush_iraq_dem_se_3, bush_iraq_dem_se_4, 
                    obama_iraq_dem_se_5, obama_iraq_dem_se_6,
                    iraq_bush_cces_dem_se_7, iraq_bush_cces_dem_se_8, 
                    iraq_obama_cces_dem_se_9, iraq_obama_cces_dem_se_10),
          covariate.labels = c("Black", "Male",
                               "Some College", "College Grad", "Post-Grad",
                               "Age: 30-49", "Age: 50-64", "Age: 65+",
                               "Income: $30K-$50K",
                               "Income: $50K-$75K",
                               "Income: $75K and above",
                               rep(NA, 61), "Constant"),
          omit = c("poll"),
          omit.labels = c("Poll FE"),
          column.labels = c("Bush (Casualties)", "Bush (Alienation)",
                            "Obama (Alienation)", "Bush (CCES)" ,"Obama (CCES)"),
          column.separate = c(2,2,2,2,2),
          add.lines = list(c("Polls", 18, 17, 39, 35,
                             11,11,2,2,3,3)),
          omit.stat = c("adj.rsq","f", "ser"),
          out = "table_b1_democrats_full_regression_new.tex")


####################################################
# 2e) Generate Figure 11 in the Appendix
####################################################
coef_black<-c(coef(bush_iraq_casualty_dem_reg_1)[2],
              coef(bush_iraq_casualty_dem_reg_2)[2],
              coef(bush_iraq_dem_reg_3)[2],
              coef(bush_iraq_dem_reg_4)[2],  
              coef(obama_iraq_dem_reg_5)[2], 
              coef(obama_iraq_dem_reg_6)[2], 
              coef(iraq_bush_cces_dem_reg_7)[2],
              coef(iraq_bush_cces_dem_reg_8)[2], 
              coef(iraq_obama_cces_dem_reg_9)[2],
              coef(iraq_obama_cces_dem_reg_10)[2])

dem_se_black<-c(bush_iraq_casualty_dem_se_1[2],
            bush_iraq_casualty_dem_se_2[2], 
            bush_iraq_dem_se_3[2],
            bush_iraq_dem_se_4[2], 
            obama_iraq_dem_se_5[2],
            obama_iraq_dem_se_6[2], 
            iraq_bush_cces_dem_se_7[2],
            iraq_bush_cces_dem_se_8[2], 
            iraq_obama_cces_dem_se_9[2],
            iraq_obama_cces_dem_se_10[2])

poll_names<-c("ABC\n(Bush, Iraq)", "ABC\n(Bush, Iraq)",
              "Asstd.\n(Bush, Iraq)", "Asstd.\n(Bush, Iraq)", 
              "Asstd.\n(Obama, Iraq)", "Asstd.\n(Obama, Iraq)",
              "CCES\n(Bush, Iraq)", "CCES\n(Bush, Iraq)",
              "CCES\n(Obama, Iraq)","CCES\n(Obama, Iraq)")

poll_type<-c("Bivariate", "Full Regression",
             "Bivariate", "Full Regression",
             "Bivariate", "Full Regression",
             "Bivariate", "Full Regression",
             "Bivariate", "Full Regression")

dat<-data.frame(names = poll_names,
                type = poll_type,
                coef = coef_black,
                se = dem_se_black)

dat$ub<-dat$coef+1.96*dat$se
dat$lb<-dat$coef-1.96*dat$se

dat$ub90<-dat$coef+1.645*dat$se
dat$lb90<-dat$coef-1.645*dat$se

dat$names=factor(dat$names, levels = c("ABC\n(Bush, Iraq)", 
                                       "Asstd.\n(Bush, Iraq)", 
                                       "Asstd.\n(Obama, Iraq)", 
                                       "CCES\n(Bush, Iraq)", "CCES\n(Obama, Iraq)"))

pdf(file="black_coef_dem_new.pdf", width = 10, height = 5)
ggplot(dat, aes(names, coef, group = type, color = type))+
  geom_errorbar(ymin = dat$lb, ymax = dat$ub, width = 0.2, 
                position=position_dodge(width=0.5))+
  geom_point(position=position_dodge(width=0.5))+
  scale_y_continuous(name = "Coefficient on 'Black'",limits = c(-0.4, 0.2),
                     breaks = c(-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2))+
  geom_hline(yintercept = 0)+scale_x_discrete(name = "")+
  scale_color_manual(name = "Specification", values = c("dodgerblue", "blue4"))+theme_bw()
dev.off()

####################################################
# 3) Graph of different polls by year
####################################################
bush_iraq_graph_dat<-data.frame(matrix(rep(NA, 12*2000), nrow = 2000, ncol = 12))
obama_iraq_graph_dat<-data.frame(matrix(rep(NA, 6*2000), nrow = 2000, ncol = 6))
bush_iraq_abc_graph_dat<-data.frame(matrix(rep(NA, 12*2000), nrow = 2000, ncol = 12))
bush_iraq_cces_graph_dat<-data.frame(matrix(rep(NA, 4*2000), nrow = 2000, ncol = 4))
obama_iraq_cces_graph_dat<-data.frame(matrix(rep(NA, 6*2000), nrow = 2000, ncol = 6))

set.seed(02138)

# Bootstrap standard errors for each poll/year/group combination
# Repeated sampling within each poll

# NOTE: Takes ~10-15 minutes to run this bootstrap 2000x

for(i in 1:2000){
  ############### Bush Asstd

  bush_iq_asstd_polls<-unique(Bush_data_IQ_BW$poll)[!is.na(unique(Bush_data_IQ_BW$poll))]
  bush_iq_asstd_index<-c()

  for(j in 1:length(bush_iq_asstd_polls)){
    b_w_sample<-c(sample(which(Bush_data_IQ_BW$poll==bush_iq_asstd_polls[j]&Bush_data_IQ_BW$black==0),
                         length(which(Bush_data_IQ_BW$poll==bush_iq_asstd_polls[j]&Bush_data_IQ_BW$black==0)),
                         replace = T),
                  sample(which(Bush_data_IQ_BW$poll==bush_iq_asstd_polls[j]&Bush_data_IQ_BW$black==1),
                         length(which(Bush_data_IQ_BW$poll==bush_iq_asstd_polls[j]&Bush_data_IQ_BW$black==1)),
                         replace = T))

    bush_iq_asstd_index<-c(bush_iq_asstd_index,
                           b_w_sample)
  }

  bush_iraq_sample<-Bush_data_IQ_BW[bush_iq_asstd_index,]

  bush_iraq_graph_dat[i,]<-aggregate(bush_iraq_sample$outcome,
                                   by = list(bush_iraq_sample$year,bush_iraq_sample$black),
                                   function(x) mean(x, na.rm = T))$x


  ################ Obama Asstd.
  obama_iq_asstd_polls<-unique(Obama_data_IQ_BW$poll)[!is.na(unique(Obama_data_IQ_BW$poll))]
  obama_iq_asstd_index<-c()

  for(k in 1:length(obama_iq_asstd_polls)){
    b_w_sample<-c(sample(which(Obama_data_IQ_BW$poll==obama_iq_asstd_polls[k]&Obama_data_IQ_BW$black==0),
                       length(which(Obama_data_IQ_BW$poll==obama_iq_asstd_polls[k]&Obama_data_IQ_BW$black==0)),
                       replace = T),
    sample(which(Obama_data_IQ_BW$poll==obama_iq_asstd_polls[k]&Obama_data_IQ_BW$black==1),
           length(which(Obama_data_IQ_BW$poll==obama_iq_asstd_polls[k]&Obama_data_IQ_BW$black==1)),
           replace = T))

    obama_iq_asstd_index<-c(obama_iq_asstd_index,
                            b_w_sample)
  }

  obama_iraq_sample<-Obama_data_IQ_BW[obama_iq_asstd_index,]

  obama_iraq_graph_dat[i,]<-aggregate(obama_iraq_sample$outcome,
                                     by = list(obama_iraq_sample$year,obama_iraq_sample$black),
                                     function(x) mean(x, na.rm = T))$x

  ######### Bush ABC
  bush_iq_abc_polls<-unique(abc_iraq_combine$poll_code)[!is.na(unique(abc_iraq_combine$poll_code))]
  bush_iq_abc_index<-c()

  for(j in 1:length(bush_iq_abc_polls)){
    b_w_sample<-c(sample(which(abc_iraq_combine$poll_code==bush_iq_abc_polls[j]&abc_iraq_combine$black==0),
                         length(which(abc_iraq_combine$poll_code==bush_iq_abc_polls[j]&abc_iraq_combine$black==0)),
                         replace = T),
                  sample(which(abc_iraq_combine$poll_code==bush_iq_abc_polls[j]&abc_iraq_combine$black==1),
                         length(which(abc_iraq_combine$poll_code==bush_iq_abc_polls[j]&abc_iraq_combine$black==1)),
                         replace = T))

    bush_iq_abc_index<-c(bush_iq_abc_index,
                           b_w_sample)
  }

  bush_iraq_abc_sample<-abc_iraq_combine[bush_iq_abc_index,]

  bush_iraq_abc_graph_dat[i,]<-aggregate(bush_iraq_abc_sample$support_war,
                                     by = list(bush_iraq_abc_sample$year,bush_iraq_abc_sample$black),
                                     function(x) mean(x, na.rm = T))$x

  print(i)

  #### Bush CCES
  bush_iraq_cces_index<-c(sample(which(cces_bush$year==2006&cces_bush$black==0),
                                 length(which(cces_bush$year==2006&cces_bush$black==0)), replace = T),
                          sample(which(cces_bush$year==2006&cces_bush$black==1),
                                 length(which(cces_bush$year==2006&cces_bush$black==1)), replace = T),
                          sample(which(cces_bush$year==2008&cces_bush$black==0),
                                 length(which(cces_bush$year==2008&cces_bush$black==0)), replace = T),
                          sample(which(cces_bush$year==2008&cces_bush$black==1),
                                 length(which(cces_bush$year==2008&cces_bush$black==1)), replace = T))
  bush_iraq_cces_sample<-cces_bush[bush_iraq_cces_index,]

  bush_iraq_cces_graph_dat[i,]<-aggregate(bush_iraq_cces_sample$iraq_support,
                                          by = list(bush_iraq_cces_sample$year,
                                                    bush_iraq_cces_sample$black),
                                          function(x) mean(x, na.rm = T))$x

  #### Obama CCES
  obama_iraq_cces_index<-c(sample(which(cces_obama$year==2010&cces_obama$black==0),
                                 length(which(cces_obama$year==2010&cces_obama$black==0)), replace = T),
                          sample(which(cces_obama$year==2010&cces_obama$black==1),
                                 length(which(cces_obama$year==2010&cces_obama$black==1)), replace = T),
                          sample(which(cces_obama$year==2012&cces_obama$black==0),
                                 length(which(cces_obama$year==2012&cces_obama$black==0)), replace = T),
                          sample(which(cces_obama$year==2012&cces_obama$black==1),
                                 length(which(cces_obama$year==2012&cces_obama$black==1)), replace = T),
                          sample(which(cces_obama$year==2014&cces_obama$black==0),
                                 length(which(cces_obama$year==2014&cces_obama$black==0)), replace = T),
                          sample(which(cces_obama$year==2014&cces_obama$black==1),
                                 length(which(cces_obama$year==2014&cces_obama$black==1)), replace = T))

  obama_iraq_cces_sample<-cces_obama[obama_iraq_cces_index,]

  obama_iraq_cces_graph_dat[i,]<-aggregate(obama_iraq_cces_sample$iraq_support,
                                          by = list(obama_iraq_cces_sample$year,
                                                    obama_iraq_cces_sample$black),
                                          function(x) mean(x, na.rm = T))$x
}

####################################################
# 3a) Generate Figure 2 (polling averages by year)
####################################################
bush_iq_year_means<-colMeans(bush_iraq_graph_dat)
bush_iq_year_ub<-apply(bush_iraq_graph_dat, 2, function(x) {quantile(x, 0.975)})
bush_iq_year_lb<-apply(bush_iraq_graph_dat, 2, function(x) {quantile(x, 0.025)})

bush_iraq_asstd_graph_dat<-data.frame(means = bush_iq_year_means,
           ub = bush_iq_year_ub,
           lb = bush_iq_year_lb,
           year = c(2002:2007, 2002:2007),
           group = c(rep("white", 6), rep("black", 6)),
           poll = rep("Bush Asstd.", 12))

bush_abc_year_means<-colMeans(bush_iraq_abc_graph_dat)
bush_abc_year_ub<-apply(bush_iraq_abc_graph_dat, 2, function(x) {quantile(x, 0.975)})
bush_abc_year_lb<-apply(bush_iraq_abc_graph_dat, 2, function(x) {quantile(x, 0.025)})

bush_abc_graph_dat<-data.frame(means = bush_abc_year_means,
                                      ub = bush_abc_year_ub,
                                      lb = bush_abc_year_lb,
                                      year = c(2003:2008, 2003:2008),
                                      group = c(rep("white", 6), rep("black", 6)),
                                      poll = rep("Bush ABC", 12))


obama_iq_year_means<-colMeans(obama_iraq_graph_dat)
obama_iq_year_ub<-apply(obama_iraq_graph_dat, 2, function(x) {quantile(x, 0.975)})
obama_iq_year_lb<-apply(obama_iraq_graph_dat, 2, function(x) {quantile(x, 0.025)})

obama_iraq_asstd_graph_dat<-data.frame(means = obama_iq_year_means,
                                       ub = obama_iq_year_ub,
                                       lb = obama_iq_year_lb,
                                       year = c(2009:2011, 2009:2011),
                                       group = c(rep("white", 3), rep("black", 3)),
                                       poll = rep("Obama Asstd.", 6))
## CCES
bush_cces_year_means<-colMeans(bush_iraq_cces_graph_dat)
bush_cces_year_ub<-apply(bush_iraq_cces_graph_dat, 2, function(x) {quantile(x, 0.975)})
bush_cces_year_lb<-apply(bush_iraq_cces_graph_dat, 2, function(x) {quantile(x, 0.025)})

bush_iraq_cces_graph_dat<-data.frame(means = bush_cces_year_means,
                                      ub = bush_cces_year_ub,
                                      lb = bush_cces_year_lb,
                                      year = c(2006, 2008,2006, 2008),
                                      group = c(rep("white", 2), rep("black", 2)),
                                      poll = rep("Bush CCES", 4))

obama_cces_year_means<-colMeans(obama_iraq_cces_graph_dat)
obama_cces_year_ub<-apply(obama_iraq_cces_graph_dat, 2, function(x) {quantile(x, 0.975)})
obama_cces_year_lb<-apply(obama_iraq_cces_graph_dat, 2, function(x) {quantile(x, 0.025)})

obama_iraq_cces_graph_dat<-data.frame(means = obama_cces_year_means,
                                      ub = obama_cces_year_ub,
                                      lb = obama_cces_year_lb,
                                      year = c(seq(2010, 2014, 2),seq(2010, 2014, 2)),
                                      group = c(rep("white", 3), rep("black", 3)),
                                      poll = rep("Obama CCES", 6))

combined_cces_graph<-rbind(bush_iraq_cces_graph_dat,
                           obama_iraq_cces_graph_dat)
#####################

combined_polls<-rbind(bush_iraq_asstd_graph_dat, obama_iraq_asstd_graph_dat)

library(cowplot)

asstd_graph<-ggplot(combined_polls, aes(year, means, group = group, color = factor(group)))+
  geom_errorbar(ymin = combined_polls$lb, ymax = combined_polls$ub, width = 0.2)+
  #geom_point()+
  geom_line()+
  scale_y_continuous("% Supporting Iraq War", breaks = seq(0, 1, 0.25),
                     labels = seq(0, 100, 25),
                     limits = c(0, 1))+
  scale_x_continuous("", breaks = seq(2002, 2014, 2),
                     limits = c(2001.5, 2014.5))+
  geom_vline(xintercept = 2008.5)+
  scale_color_manual("Group", labels = c("Black", "White"),
                     values = c("black", "grey55"))+theme_bw()+theme_set(theme_bw(base_size = 10, base_family = "EB Garamond" ))+ggtitle("Assorted Polls (Alienation)")

common_legend<-get_legend(asstd_graph)
asstd_graph<-asstd_graph+theme(legend.position="none")

abc_graph<-ggplot(bush_abc_graph_dat, aes(year, means, group = group, color = factor(group)))+
  geom_errorbar(ymin = bush_abc_graph_dat$lb, ymax = bush_abc_graph_dat$ub, width = 0.2)+
  #geom_point()+
  geom_line()+
  scale_y_continuous("% Supporting Iraq War", breaks = seq(0, 1, 0.25),
                     labels = seq(0, 100, 25),
                     limits = c(0, 1))+
  scale_x_continuous("", breaks = seq(2002, 2014, 2),
                     limits = c(2001.5, 2014.5))+
  geom_vline(xintercept = 2008.5)+ggtitle("ABC Polls (Casualties)")+
  scale_color_manual("Group", values = c("black", "grey55"))+theme_bw()+theme_set(theme_bw(base_size = 10, base_family = "EB Garamond" ))+
    theme(legend.position="none")

cces_graph<-ggplot(combined_cces_graph, aes(year, means, group = group, color = factor(group)))+
  geom_errorbar(ymin = combined_cces_graph$lb, ymax = combined_cces_graph$ub, width = 0.2, show.legend = FALSE)+
  #geom_point()+
  geom_line()+
  scale_y_continuous("% Supporting Iraq War", breaks = seq(0, 1, 0.25),
                     labels = seq(0, 100, 25),
                     limits = c(0, 1))+
  scale_x_continuous("", breaks = seq(2002, 2014, 2),
                     limits = c(2001.5, 2014.5))+
  geom_vline(xintercept = 2008.5)+ggtitle("CCES Polls")+
  scale_color_manual("Group", values = c("black", "grey55"))+
  theme_bw()+theme_set(theme_bw(base_size = 10, base_family = "EB Garamond"))+theme(legend.position="none")

library(gridExtra)

pdf("all_datasets_gap_new_garamond.pdf", width = 6, height = 5)
grid.arrange(abc_graph, asstd_graph, cces_graph,
             common_legend,
             ncol=2, widths=c(2, 2))
dev.off()


setEPS()                                             
postscript("all_datasets_gap_new_garamond.eps", width = 6, height = 5) 
grid.arrange(abc_graph, asstd_graph, cces_graph,
             common_legend,
             ncol=2, widths=c(2, 2))
dev.off()

###################################################################
#
# 4) Risks of War - Casualties & Military Service
#
###################################################################

####################################################################################
# 4a) Regressions for Figure 4 & Table A4, Models 1-3 (Casualties - Direct Effects)
####################################################################################

# Model 1 - correlates of casualty sensitivity
bush_iq_cas_reg_1<-lm(casual_sens ~ black+gender+party_id+education+age+income+
                               factor(poll), 
                             data = abc_iraq_combine)

bush_iq_cas_se_1<-sqrt(diag(vcovHC(bush_iq_cas_reg_1, type = "HC1")))
coeftest(bush_iq_cas_reg_1, vcovHC(bush_iq_cas_reg_1, type = "HC1"))

# Model 2 - support for war, black sample only
bush_iq_cas_reg_2<-lm(support_war ~ casual_sens+gender+party_id+education+age+income+
                        factor(poll), 
                      data = abc_iraq_combine[which(abc_iraq_combine$black==1),])

bush_iq_cas_se_2<-sqrt(diag(vcovHC(bush_iq_cas_reg_2, type = "HC1")))
coeftest(bush_iq_cas_reg_2, vcovHC(bush_iq_cas_reg_2, type = "HC1"))

# values
coef(bush_iq_cas_reg_1)[2]+1.96*bush_iq_cas_se_1[2]
coef(bush_iq_cas_reg_1)[2]-1.96*bush_iq_cas_se_1[2]

# Model 3 - support for war, white sample only
bush_iq_cas_reg_3<-lm(support_war ~ casual_sens+gender+party_id+education+age+income+
                        factor(poll), 
                      data = abc_iraq_combine[which(abc_iraq_combine$black==0),])

bush_iq_cas_se_3<-sqrt(diag(vcovHC(bush_iq_cas_reg_3, type = "HC1")))
coeftest(bush_iq_cas_reg_3, vcovHC(bush_iq_cas_reg_3, type = "HC1"))

# Mediation Analysis - Casualty Sensitivity
abc_med <- abc_iraq_combine[-which(is.na(abc_iraq_combine$support_war)),]

set.seed(02138)

med.fit.bush.cas <- lm(casual_sens ~ black+gender+income+education+party_id+age+
                factor(poll), data=abc_med)
out.fit.bush.cas <- lm(support_war ~ black + casual_sens+gender+income+education+party_id+age+
                factor(poll), data=abc_med)
med.out.bush.cas <- mediate(med.fit.bush.cas, out.fit.bush.cas, treat = "black", mediator = "casual_sens",
                   robustSE = TRUE, sims = 250)
summary(med.out.bush.cas)
plot(med.out.bush.cas, xlim = c(-0.25, 0.2), family = "EB Garamond",)

###############################################################################
# 4b) Regressions for Figure 6 & Table A4, Models 4-6 (Military Service, Bush)
###############################################################################

# Model 4 - Correlates of military service
bush_iq_serv_reg_1<-lm(serving ~ black+gender+party_id+education+age+income+
                         factor(poll), data = cces_bush)

coeftest(bush_iq_serv_reg_1, vcovHC(bush_iq_serv_reg_1, 
                                    type = "HC1"))

bush_iq_serv_se_1<-sqrt(diag(vcovHC(bush_iq_serv_reg_1, 
                                 type = "HC1")))

# Model 5 - Correlates of support for war, black sample
bush_iq_serv_reg_2<-lm(iraq_support ~ serving+gender+party_id+education+age+income+
                         factor(poll), data = cces_bush[which(cces_bush$black_aa==1),])

coeftest(bush_iq_serv_reg_2, vcovHC(bush_iq_serv_reg_2, 
                                type = "HC1"))

bush_iq_serv_se_2<-sqrt(diag(vcovHC(bush_iq_serv_reg_2, 
                                 type = "HC1")))
# Model 6 - Correlates of support for war, white sample
bush_iq_serv_reg_3<-lm(iraq_support ~ serving+gender+party_id+education+age+income+
                         factor(poll), data = cces_bush[which(cces_bush$black_aa==0),])

coeftest(bush_iq_serv_reg_3, vcovHC(bush_iq_serv_reg_3, 
                                type = "HC1"))

bush_iq_serv_se_3<-sqrt(diag(vcovHC(bush_iq_serv_reg_3, 
                                 type = "HC1")))


# Mediation Analysis - Military Service, Bush
set.seed(02138)

med.fit.bush.serv <- lm(serving ~ black+gender+income+education+party_id+age+factor(poll), data = cces_bush[-which(is.na(cces_bush$iraq_support)),])
out.fit.bush.serv <- lm(iraq_support ~ black + serving+gender+income+education+party_id+age+
                          factor(poll), data = cces_bush[-which(is.na(cces_bush$iraq_support)),])
med.out.bush.serv <- mediate(med.fit.bush.serv, out.fit.bush.serv, treat = "black", mediator = "serving",
                   robustSE = TRUE, sims = 250)
summary(med.out.bush.serv)
plot(med.out.bush.serv, xlim = c(-0.25, 0.2),  family = "EB Garamond")

###################################################################
# 4c) Regressions for Figure 6 & Table A4, Models 7-9 (Military Service, Obama)
###################################################################

# Model 7 - Correlates of military service
obama_iq_serv_reg_1<-lm(serving ~ black+gender+party_id+education+age+income+
                          factor(poll), data = cces_obama)

coeftest(obama_iq_serv_reg_1, vcovHC(obama_iq_serv_reg_1, 
                                     type = "HC1"))

obama_iq_serv_se_1<-sqrt(diag(vcovHC(obama_iq_serv_reg_1, 
                                       type = "HC1")))

# Model 8 - Correlates of support for war, black sample
obama_iq_serv_reg_2<-lm(iraq_support ~ serving+gender+party_id+education+age+income+
                          factor(poll), 
                        data = cces_obama[which(cces_obama$black_aa==1),])

coeftest(obama_iq_serv_reg_2, vcovHC(obama_iq_serv_reg_2, 
                                       type = "HC1"))

obama_iq_serv_se_2<-sqrt(diag(vcovHC(obama_iq_serv_reg_2, 
                                       type = "HC1")))

# Model 9 - Correlates of support for war, white sample
obama_iq_serv_reg_3<-lm(iraq_support ~ serving+gender+party_id+education+age+income+
                          factor(poll), 
                          data = cces_obama[which(cces_obama$black_aa==0),])

coeftest(obama_iq_serv_reg_3, vcovHC(obama_iq_serv_reg_3, 
                                       type = "HC1"))

obama_iq_serv_se_3<-sqrt(diag(vcovHC(obama_iq_serv_reg_3, 
                                       type = "HC1")))

# values
coef(bush_iq_serv_reg_1)[2]+1.96*bush_iq_serv_se_1[2]
coef(bush_iq_serv_reg_1)[2]-1.96*bush_iq_serv_se_1[2]

coef(obama_iq_serv_reg_1)[2]+1.96*obama_iq_serv_se_1[2]
coef(obama_iq_serv_reg_1)[2]-1.96*obama_iq_serv_se_1[2]

# Mediation Analysis - Military Service, Obama
set.seed(02138)

med.fit.obama.serv <- lm(serving ~ black + 
                gender + party_id + education +
                age + income + factor(poll), data = cces_obama[-which(is.na(cces_obama$iraq_support)),])
out.fit.obama.serv <- lm(iraq_support ~ black + serving + 
                gender + party_id + education +
                age + income + factor(poll), data = cces_obama[-which(is.na(cces_obama$iraq_support)),])
med.out.obama.serv <- mediate(med.fit.obama.serv, out.fit.obama.serv, 
                              treat = "black", mediator = "serving",
                   robustSE = TRUE, sims = 250)
summary(med.out.obama.serv)
plot(med.out.obama.serv, xlim = c(-0.25, 0.2),  family = "EB Garamond")

###################################################################
# 5d) Figure 6, Mediation Plot for Costs of War
###################################################################
pdf("casualty_mediation_plots_new_garamond.pdf", width = 4, height = 5)
par(mar=c(2, 4, 4, 2),
        mfcol = c(3,1))
plot(med.out.bush.cas, xlim = c(-0.25, 0.2),  family = "EB Garamond")
title(main = list("Casualty Sensitivity (Bush)", cex = 1.5,
                  col = "black", font = 1), line = 1, adj = 0,  family = "EB Garamond")
plot(med.out.bush.serv, xlim = c(-0.25, 0.2),  family = "EB Garamond")
title(main = list("Military Family (Bush)", cex = 1.5,
                  col = "black", font = 1), line = 1, adj = 0,  family = "EB Garamond")
plot(med.out.obama.serv, xlim = c(-0.25, 0.2),  family = "EB Garamond")
title(main = list("Military Family (Obama)", cex = 1.5,
                  col = "black", font = 1), line = 1, adj = 0,  family = "EB Garamond")
dev.off()

setEPS()                                             
postscript("casualty_mediation_plots_new_garamond.eps", width = 4, height = 5)
par(mar=c(2, 4, 4, 2),
    mfcol = c(3,1))
plot(med.out.bush.cas, xlim = c(-0.25, 0.2),  family = "EB Garamond")
title(main = list("Casualty Sensitivity (Bush)", cex = 1.5,
                  col = "black", font = 1), line = 1, adj = 0,  family = "EB Garamond")
plot(med.out.bush.serv, xlim = c(-0.25, 0.2),  family = "EB Garamond")
title(main = list("Military Family (Bush)", cex = 1.5,
                  col = "black", font = 1), line = 1, adj = 0,  family = "EB Garamond")
plot(med.out.obama.serv, xlim = c(-0.25, 0.2),  family = "EB Garamond")
title(main = list("Military Family (Obama)", cex = 1.5,
                  col = "black", font = 1), line = 1, adj = 0,  family = "EB Garamond")
dev.off()

###################################################################
#5e) Generates Table 3, Models 1-9 (Costs of War)
###################################################################
stargazer(bush_iq_cas_reg_1,
          bush_iq_cas_reg_2,
          bush_iq_cas_reg_3,
          bush_iq_serv_reg_1,
          bush_iq_serv_reg_2,
          bush_iq_serv_reg_3,
          obama_iq_serv_reg_1,
          obama_iq_serv_reg_2,
          obama_iq_serv_reg_3,
          se = list(bush_iq_cas_se_1,
                    bush_iq_cas_se_2,
                    bush_iq_cas_se_3,
                    bush_iq_serv_se_1,
                    bush_iq_serv_se_2,
                    bush_iq_serv_se_3,
                    obama_iq_serv_se_1,
                    obama_iq_serv_se_2,
                    obama_iq_serv_se_3),
          column.labels = rep(c("Combined", "Black", "White"), 3),
          covariate.labels = c("Black", "Casualty Sens.", "Military Family", "Male",
                               "Independent", "Republican",
                               "Some College", "College Grad", "Post-Grad",
                               "Age: 30-49", "Age: 50-64", "Age: 65+",
                               "Income: $30K-$50K",
                               "Income: $50K-$75K",
                               "Income: $75K and above",
                               rep(NA, 19), "Constant"),
          omit = c("poll"),
          omit.labels = c("Poll FE"),
          add.lines = list(c("Polls", 17, 17, 17, 2,2,2,3,3,3)),
          omit.stat = c("adj.rsq","f", "ser"),
          out = "table_A4_cas_regression_new.tex")

###################################################################
#
# 6) Costs of War - Spending
#
###################################################################

###################################################################
# 6a) Regressions for Table 4, Models 1-3 (Spending, Bush)
###################################################################

# Model 1 - correlates of support for military spending, Bush
bush_iq_spend_reg_1<-lm(dom_over_def ~ black+gender+party_id+education+age+income+
                          factor(poll), data = cces_bush)

coeftest(bush_iq_spend_reg_1, vcovHC(bush_iq_spend_reg_1, 
                                    type = "HC1"))

bush_iq_spend_se_1<-sqrt(diag(vcovHC(bush_iq_spend_reg_1, 
                                    type = "HC1")))

# Model 2 - correlates of support for war, black sample
bush_iq_spend_reg_2<-lm(iraq_support ~ dom_over_def+gender+party_id+education+age+income+
                          factor(poll), data = cces_bush[which(cces_bush$black_aa==1),])

coeftest(bush_iq_spend_reg_2, vcovHC(bush_iq_spend_reg_2, 
                                    type = "HC1"))

bush_iq_spend_se_2<-sqrt(diag(vcovHC(bush_iq_spend_reg_2, 
                                    type = "HC1")))

# Model 3 - correlates of support for war, white sample
bush_iq_spend_reg_3<-lm(iraq_support ~ dom_over_def+gender+party_id+education+age+income+
                          factor(poll), data = cces_bush[which(cces_bush$black_aa==0),])

coeftest(bush_iq_spend_reg_3, vcovHC(bush_iq_spend_reg_3, 
                                    type = "HC1"))

bush_iq_spend_se_3<-sqrt(diag(vcovHC(bush_iq_spend_reg_3, 
                                    type = "HC1")))


# Mediation Analysis - Spending, Bush
set.seed(02138)

med.fit.bush.spend <- lm(dom_over_def ~ black + 
                          gender + party_id + education +
                          age + income + factor(poll), data = cces_bush[-which(is.na(cces_bush$iraq_support)),])
out.fit.bush.spend <- lm(iraq_support ~ black + dom_over_def + 
                          gender + party_id + education +
                          age + income + factor(poll), data = cces_bush[-which(is.na(cces_bush$iraq_support)),])
med.out.bush.spend <- mediate(med.fit.bush.spend, out.fit.bush.spend, treat = "black", mediator = "dom_over_def",
                             robustSE = TRUE, sims = 250)
summary(med.out.bush.spend)
plot(med.out.bush.spend, xlim = c(-0.25, 0.2),  family = "EB Garamond")

###################################################################
# 6b) Regressions for Table 5, Models 4-6 (Spending, Obama)
###################################################################

# Model 4 - correlates of support for military spending, Obama
obama_iq_spend_reg_1<-lm(dom_over_def ~ black+gender+party_id+education+age+income+
                           factor(poll), data = cces_obama)

coeftest(obama_iq_spend_reg_1, vcovHC(obama_iq_spend_reg_1, 
                                     type = "HC1"))

obama_iq_spend_se_1<-sqrt(diag(vcovHC(obama_iq_spend_reg_1, 
                                     type = "HC1")))

# Model 5 - correlates of support for war, black sample
obama_iq_spend_reg_2<-lm(iraq_support ~ dom_over_def+gender+party_id+education+age+income+
                           factor(poll), data = cces_obama[which(cces_obama$black_aa==1),])

coeftest(obama_iq_spend_reg_2, vcovHC(obama_iq_spend_reg_2, 
                                     type = "HC1"))

obama_iq_spend_se_2<-sqrt(diag(vcovHC(obama_iq_spend_reg_2, 
                                     type = "HC1")))

# Model 6 - correlates of support for war, white sample
obama_iq_spend_reg_3<-lm(iraq_support ~ dom_over_def+gender+party_id+education+age+income+
                           factor(poll), 
                        data = cces_obama[which(cces_obama$black_aa==0),])

coeftest(obama_iq_spend_reg_3, vcovHC(obama_iq_spend_reg_3, 
                                     type = "HC1"))

obama_iq_spend_se_3<-sqrt(diag(vcovHC(obama_iq_spend_reg_3, 
                                     type = "HC1")))

# values
coef(bush_iq_spend_reg_1)[2]+1.96*bush_iq_spend_se_1[2]
coef(bush_iq_spend_reg_1)[2]-1.96*bush_iq_spend_se_1[2]

coef(obama_iq_spend_reg_1)[2]+1.96*obama_iq_spend_se_1[2]
coef(obama_iq_spend_reg_1)[2]-1.96*obama_iq_spend_se_1[2]

# Mediation Analysis - Spending, Obama
set.seed(02138)

med.fit.obama.spend <- lm(dom_over_def ~ black + 
                           gender + party_id + education +
                           age + income + factor(poll), data = cces_obama[-which(is.na(cces_obama$iraq_support)),])
out.fit.obama.spend <- lm(iraq_support ~ black + dom_over_def + 
                           gender + party_id + education +
                           age + income + factor(poll), data = cces_obama[-which(is.na(cces_obama$iraq_support)),])
med.out.obama.spend <- mediate(med.fit.obama.spend, out.fit.obama.spend, treat = "black", mediator = "dom_over_def",
                              robustSE = TRUE, sims = 250)
summary(med.out.obama.spend)
plot(med.out.obama.spend, xlim = c(-0.25, 0.2),  family = "EB Garamond")

###################################################################
# 6c) Generates Table 4, Models 1-6 (Spending & War)
###################################################################
stargazer(bush_iq_spend_reg_1,
          bush_iq_spend_reg_2,
          bush_iq_spend_reg_3,
          obama_iq_spend_reg_1,
          obama_iq_spend_reg_2,
          obama_iq_spend_reg_3,
          se = list(bush_iq_spend_se_1,
                    bush_iq_spend_se_2,
                    bush_iq_spend_se_3,
                    obama_iq_spend_se_1,
                    obama_iq_spend_se_2,
                    obama_iq_spend_se_3),
          covariate.labels = c("Black", "Dom. Spending",
                               "Male",
                               "Independent", "Republican",
                               "Some College", "College Grad", "Post-Grad",
                               "Age: 30-49", "Age: 50-64", "Age: 65+",
                               "Income: $30K-$50K",
                               "Income: $50K-$75K",
                               "Income: $75K and above",
                               rep(NA, 3), "Constant"),
          column.labels = rep(c("Combined", "Black", "White"),2),
          omit = c("poll"),
          omit.labels = c("Poll FE"),
          add.lines = list(c("Polls", 2,2,2,3,3,3)),
          omit.stat = c("adj.rsq","f", "ser"),
          out = "table_A5_spending_regression_new.tex")

###################################################################
# 6d) Figure 8, Mediation Plot for Spending & War
###################################################################
pdf("spending_mediation_plots_new_garamond.pdf", width = 4, height = 4)
par(mar=c(2, 4, 4, 2),
    mfcol = c(2, 1))
plot(med.out.bush.spend, xlim = c(-0.2, 0.2), cex.axis = .8,  family = "EB Garamond")
title(main = list("Spending (Bush)", cex = 1,
                  col = "black", font = 1), line = 1, adj = 0,  family = "EB Garamond")
plot(med.out.obama.spend, xlim = c(-0.2, 0.2), cex.axis = .8,  family = "EB Garamond")
title(main = list("Spending (Obama)", cex = 1,
                  col = "black", font = 1), line = 1, adj = 0,  family = "EB Garamond")
dev.off()

setEPS()                                             
postscript("spending_mediation_plots_new_garamond.eps", width = 4, height = 4)
par(mar=c(2, 4, 4, 2),
    mfcol = c(2, 1))
plot(med.out.bush.spend, xlim = c(-0.2, 0.2), cex.axis = .8,  family = "EB Garamond")
title(main = list("Spending (Bush)", cex = 1,
                  col = "black", font = 1), line = 1, adj = 0,  family = "EB Garamond")
plot(med.out.obama.spend, xlim = c(-0.2, 0.2), cex.axis = .8,  family = "EB Garamond")
title(main = list("Spending (Obama)", cex = 1,
                  col = "black", font = 1), line = 1, adj = 0,  family = "EB Garamond")
dev.off()

###################################################################
#
# 7) Prerequisites for War - Alienation
#
###################################################################

#############################################################################
# 7a) Regressions for Figure 9 & Table A6, Models 1-3 (Alienation, Bush)
#############################################################################

# Model 1 - correlates of alienation, Bush
Bush1alien_IQ <- lm(alienation_binary ~ 
                      black+gender+party_id+education+age+income+factor(poll), data = Bush_data_IQ_BW)
Bush1alien_IQ_se<-sqrt(diag(vcovHC(Bush1alien_IQ, type = "HC1")))

# Model 2 - correlates of support for war, black sample
Bush2alien_IQ <- lm(outcome ~ alienation_binary+gender+party_id+education+age+income+
                      poll, data = Bush_data_IQ_BW[which(Bush_data_IQ_BW$black==1),])
Bush2alien_IQ_se<-sqrt(diag(vcovHC(Bush2alien_IQ, type = "HC1")))

# Model 3 - correlates of support for war, white sample
Bush3alien_IQ <- lm(outcome ~ alienation_binary+gender+party_id+education+age+income+
                      poll, data = Bush_data_IQ_BW[which(Bush_data_IQ_BW$black==0),])
Bush3alien_IQ_se<-sqrt(diag(vcovHC(Bush3alien_IQ, type = "HC1")))

coef(Bush1alien_IQ)[2]+1.96*Bush1alien_IQ_se[2]
coef(Bush1alien_IQ)[2]-1.96*Bush1alien_IQ_se[2]
  
# Mediation - Bush & Alienation
set.seed(02138)

med.fit.bush.alien <- lm(alienation_binary ~ black + 
                           gender+party_id+education+age+income+factor(poll), data = Bush_data_IQ_BW[!is.na(Bush_data_IQ_BW$outcome),])
out.fit.bush.alien <- lm(outcome ~ black + alienation_binary + 
                           gender+party_id+education+age+income+factor(poll), data = Bush_data_IQ_BW[!is.na(Bush_data_IQ_BW$outcome),])
med.out.bush.alien <- mediate(med.fit.bush.alien, out.fit.bush.alien, 
                              treat = "black", mediator = "alienation_binary",
                              robustSE = TRUE, sims = 250)
summary(med.out.bush.alien)
plot(med.out.bush.alien, xlim = c(-0.25, 0.2),  family = "EB Garamond")


##############################################################################
# 7b) Regressions for Figure 9 & Table A6, Models 4-6 (Alienation, Obama)
##############################################################################

# Model 4 - correlates of alienation, Obama
Obama1alien_IQ <- lm(alienation_binary ~ 
                       black+gender+party_id+education+age+income+factor(poll), data = Obama_data_IQ_BW)
Obama1alien_IQ_se<-sqrt(diag(vcovHC(Obama1alien_IQ, type = "HC1")))

# Model 5 - correlates of support for war, black sample
Obama2alien_IQ <- lm(outcome ~ alienation_binary +
                       +gender+party_id+education+age+income+factor(poll), data = Obama_data_IQ_BW[which(Obama_data_IQ_BW$black==1),])
Obama2alien_IQ_se<-sqrt(diag(vcovHC(Obama2alien_IQ, type = "HC1")))

# Model 6 - correlates of support for war, white sample
Obama3alien_IQ <- lm(outcome ~ alienation_binary +
                       gender+party_id+education+age+income+factor(poll), data = Obama_data_IQ_BW[which(Obama_data_IQ_BW$black==0),])
Obama3alien_IQ_se<-sqrt(diag(vcovHC(Obama3alien_IQ, type = "HC1")))

summary(Obama3alien_IQ)

#
coef(Obama1alien_IQ)[2]+1.96*Obama1alien_IQ_se[2]
coef(Obama1alien_IQ)[2]-1.96*Obama1alien_IQ_se[2]

# Mediation - Obama & Alienation
set.seed(02138)

med.fit.obama.alien <- lm(alienation_binary ~ black + 
                            gender+party_id+education+age+income+factor(poll), 
                          data = Obama_data_IQ_BW[which(!is.na(Obama_data_IQ_BW$outcome)),])
out.fit.obama.alien <- lm(outcome ~ black + alienation_binary + 
                            gender+party_id+education+age+income+factor(poll), 
                          data = Obama_data_IQ_BW[which(!is.na(Obama_data_IQ_BW$outcome)),])
med.out.obama.alien <- mediate(med.fit.obama.alien, out.fit.obama.alien, 
                               treat = "black", mediator = "alienation_binary",
                               robustSE = TRUE, sims = 250)
summary(med.out.obama.alien)
plot(med.out.obama.alien, xlim = c(-0.25, 0.2),  family = "EB Garamond")

###################################################################
# 7c) Generates Table 5, Mediation Plot for Spending & War
###################################################################
stargazer(Bush1alien_IQ,Bush2alien_IQ,Bush3alien_IQ,
          Obama1alien_IQ,Obama2alien_IQ,Obama3alien_IQ,
          se = list(Bush1alien_IQ_se,Bush2alien_IQ_se,Bush3alien_IQ_se,
                    Obama1alien_IQ_se,Obama2alien_IQ_se,Obama3alien_IQ_se),
          covariate.labels = c("Black", "Alienation",
                               "Male",
                               "Independent", "Republican",
                               "Some College", "College Grad", "Post-Grad",
                               "Age: 30-49", "Age: 50-64", "Age: 65+",
                               "Income: $30K-$50K",
                               "Income: $50K-$75K",
                               "Income: $75K and above",
                               rep(NA, 37), "Constant"),
          column.labels = rep(c("Combined", "Black", "White"),2),
          omit = c("poll"),
          omit.labels = c("Poll FE"),
          add.lines = list(c("Polls", 35,35,35,4,4,4)),
          omit.stat = c("adj.rsq","f", "ser"),
          out = "table_A6_alienation_regression6_new.tex")

###################################################################
# 7d) Figure 10, Mediation Plot for Alienation & War
###################################################################
pdf("alien_mediation_plots_new_garamond.pdf", width = 4, height = 4)
par(mar=c(2, 4, 4, 2),
    mfcol = c(2, 1))
plot(med.out.bush.alien, xlim = c(-0.2, 0.2), cex.axis = .8,  family = "EB Garamond")
title(main = list("Alienation (Bush)", cex = 1,
                  col = "black", font = 1), line = 1, adj = 0,  family = "EB Garamond")
plot(med.out.obama.alien, xlim = c(-0.2, 0.2), cex.axis = .8,  family = "EB Garamond")
title(main = list("Alienation (Obama)", cex = 1,
                  col = "black", font = 1), line = 1, adj = 0,  family = "EB Garamond")
dev.off()

setEPS()                                             
postscript("alien_mediation_plots_new_garamond.eps", width = 4, height = 4)
par(mar=c(2, 4, 4, 2),
    mfcol = c(2, 1))
plot(med.out.bush.alien, xlim = c(-0.2, 0.2), cex.axis = .8,  family = "EB Garamond")
title(main = list("Alienation (Bush)", cex = 1,
                  col = "black", font = 1), line = 1, adj = 0,  family = "EB Garamond")
plot(med.out.obama.alien, xlim = c(-0.2, 0.2), cex.axis = .8,  family = "EB Garamond")
title(main = list("Alienation (Obama)", cex = 1,
                  col = "black", font = 1), line = 1, adj = 0,  family = "EB Garamond")
dev.off()
###################################################################
# 7e) Return statistics on % alienated by poll, president
###################################################################

# % black/white alienated under Bush, Iraq War data
mean(Bush_data_IQ_BW[Bush_data_IQ_BW$black == 1 & !is.na(Bush_data_IQ_BW$black),]$alienation_binary, na.rm = T)
mean(Bush_data_IQ_BW[Bush_data_IQ_BW$black == 0 & !is.na(Bush_data_IQ_BW$black),]$alienation_binary, na.rm = T)

# % black/white alienated under Obama, Iraq War data
mean(Obama_data_IQ_BW[Obama_data_IQ_BW$black == 1 & !is.na(Obama_data_IQ_BW$black),]$alienation_binary, na.rm = T)
mean(Obama_data_IQ_BW[Obama_data_IQ_BW$black == 0 & !is.na(Obama_data_IQ_BW$black),]$alienation_binary, na.rm = T)

# % white Democrats alienated under Obama, Iraq War data
mean(Obama_data_IQ_BW[Obama_data_IQ_BW$black == 0 & !is.na(Obama_data_IQ_BW$black) & Obama_data_IQ_BW$party_id == "democrat" & !is.na(Obama_data_IQ_BW$party_id),]$alienation_binary, na.rm = T)

# All Polls, aggregates and plot
alienation_by_poll<-aggregate(Bush_data_IQ_BW$alienation_binary, 
                              by = list(Bush_data_IQ_BW$poll, 
                                        Bush_data_IQ_BW$black,
                                        Bush_data_IQ_BW$party_id), function(x) mean(x, na.rm = T))
hist(alienation_by_poll$x[which(alienation_by_poll$Group.2==1)], breaks = 20)

median(alienation_by_poll$x[which(alienation_by_poll$Group.2==1)])
mean(alienation_by_poll$x[which(alienation_by_poll$Group.2==1)])

median(alienation_by_poll$x[which(alienation_by_poll$Group.2==0&alienation_by_poll$Group.3=="democrat")])
mean(alienation_by_poll$x[which(alienation_by_poll$Group.2==0&alienation_by_poll$Group.3=="democrat")])

alienation_by_poll<-aggregate(Obama_data_IQ_BW$alienation_binary, 
                              by = list(Obama_data_IQ_BW$poll, 
                                        Obama_data_IQ_BW$black,
                                        Obama_data_IQ_BW$party_id), function(x) mean(x, na.rm = T))
hist(alienation_by_poll$x[which(alienation_by_poll$Group.2==1)], breaks = )

alienation_by_poll<-aggregate(Obama_data_IQ_BW$outcome, 
                              by = list(Obama_data_IQ_BW$black,
                                        Obama_data_IQ_BW$alienation_binary), function(x) mean(x, na.rm = T))

alienation_by_poll

###################################################################
# 7f) Alienation vs. Caring about Black People, 2005 CNN poll
###################################################################

# Of black respondents who said that Bush did not care about them, how many also said that Bush didn't care about black people?
nrow(Bush_IQ_CNN2005[Bush_IQ_CNN2005$race == "black" & !is.na(Bush_IQ_CNN2005$race) & Bush_IQ_CNN2005$alienation_binary == 1 & !is.na(Bush_IQ_CNN2005$alienation_binary) & Bush_IQ_CNN2005$black_people == 1 & !is.na(Bush_IQ_CNN2005$black_people),])/nrow(Bush_IQ_CNN2005[Bush_IQ_CNN2005$race == "black" & !is.na(Bush_IQ_CNN2005$race) & Bush_IQ_CNN2005$alienation_binary == 1 & !is.na(Bush_IQ_CNN2005$alienation_binary),])

nrow(Bush_IQ_CNN2005[Bush_IQ_CNN2005$race == "black" & !is.na(Bush_IQ_CNN2005$race) & Bush_IQ_CNN2005$alienation_binary == 1 & !is.na(Bush_IQ_CNN2005$alienation_binary) & Bush_IQ_CNN2005$black_people == 1 & !is.na(Bush_IQ_CNN2005$black_people),])
nrow(Bush_IQ_CNN2005[Bush_IQ_CNN2005$race == "black" & !is.na(Bush_IQ_CNN2005$race) & Bush_IQ_CNN2005$alienation_binary == 1 & !is.na(Bush_IQ_CNN2005$alienation_binary),])

# Of white respondents who said that Bush did not care about them, how many also said that Bush didn't care about black people?
nrow(Bush_IQ_CNN2005[Bush_IQ_CNN2005$race == "white" & !is.na(Bush_IQ_CNN2005$race) & Bush_IQ_CNN2005$alienation_binary == 1 & !is.na(Bush_IQ_CNN2005$alienation_binary) & Bush_IQ_CNN2005$black_people == 1 & !is.na(Bush_IQ_CNN2005$black_people),])/nrow(Bush_IQ_CNN2005[Bush_IQ_CNN2005$race == "white" & !is.na(Bush_IQ_CNN2005$race) & Bush_IQ_CNN2005$alienation_binary == 1 & !is.na(Bush_IQ_CNN2005$alienation_binary),])

# Of black respondents, how many said that Bush didn't care about people like them?
nrow(Bush_IQ_CNN2005[Bush_IQ_CNN2005$race == "black" & !is.na(Bush_IQ_CNN2005$race) & Bush_IQ_CNN2005$alienation_binary == 1 & !is.na(Bush_IQ_CNN2005$alienation_binary),])/nrow(Bush_IQ_CNN2005[Bush_IQ_CNN2005$race == "black" & !is.na(Bush_IQ_CNN2005$race),])

# Of white respondents, how many said that Bush didn't care about people like them?
nrow(Bush_IQ_CNN2005[Bush_IQ_CNN2005$race == "white" & !is.na(Bush_IQ_CNN2005$race) & Bush_IQ_CNN2005$alienation_binary == 1 & !is.na(Bush_IQ_CNN2005$alienation_binary),])/nrow(Bush_IQ_CNN2005[Bush_IQ_CNN2005$race == "white" & !is.na(Bush_IQ_CNN2005$race),])

####################################################
#
# 8) Plot of all explanatory variables
#
####################################################

#####################################################################################
# Output for Figures 4 & 5 (Coefficient Plot, Casualty Sensitivity & Military Family)
#####################################################################################

# Coefficients from each regression
effects_cas<-c(coef(bush_iq_cas_reg_1)[2], # bush, casualty sens
               coef(bush_iq_serv_reg_1)[2], # bush, family member serving
               coef(obama_iq_serv_reg_1)[2], # obama, family member serving
               coef(bush_iq_cas_reg_2)[2], # bush, casualty sens, black sample
               coef(bush_iq_serv_reg_2)[2], # bush, family member serving, black sample
               coef(obama_iq_serv_reg_2)[2], # obama, family member serving, black sample
               coef(bush_iq_cas_reg_3)[2], # bush, casualty sens, white sample
               coef(bush_iq_serv_reg_3)[2], # bush, family member serving, white sample
               coef(obama_iq_serv_reg_3)[2]) # obama, family member serving, white sample

# Standard errors from each regression
se_cas<-c(bush_iq_cas_se_1[2], # bush, casualty sens
               bush_iq_serv_se_1[2], # bush, family member serving
               obama_iq_serv_se_1[2], # obama, family member serving
               bush_iq_cas_se_2[2], # bush, casualty sens, black sample
               bush_iq_serv_se_2[2], # bush, family member serving, black sample
               obama_iq_serv_se_2[2], # obama, family member serving, black sample
               bush_iq_cas_se_3[2], # bush, casualty sens, white sample
               bush_iq_serv_se_3[2], # bush, family member serving, white sample
               obama_iq_serv_se_3[2]) # obama, family member serving, white sample

# Labels for each dataset used
poll_names_cas<-c(c("1A\nCas. Sensitivity\n(Bush, Iraq)"),
                  c("2A\nMil. Family\n(Bush, Iraq)"),c("2B\nMil. Family\n(Obama, Iraq)"),
                  c("1A\nCas. Sensitivity\n(Bush, Iraq)"),
                  c("2A\nMil. Family\n(Bush, Iraq)"),c("2B\nMil. Family\n(Obama, Iraq)"),
                  c("1A\nCas. Sensitivity\n(Bush, Iraq)"),
                  c("2A\nMil. Family\n(Bush, Iraq)"),c("2B\nMil. Family\n(Obama, Iraq)"))

# Category of each coefficient
effect_type_cas<-c(rep("Black (All)", 3),
                   rep("Cas. Sens. (Black)",1),rep("Mil. Family (Black)",2),
                   rep("Cas. Sens. (White)", 1),rep("Mil. Family (White)",2))

# Administration
admin_type_cas<-c(rep("Bush", 2), "Obama",
                  rep("Bush", 2), "Obama",
                  rep("Bush", 2), "Obama")

# Combine to one data frame
dat_cas<-data.frame(names = effect_type_cas,
                poll = poll_names_cas,
                admin = admin_type_cas,
                coef = effects_cas,
                se = se_cas)

# Generate upper and lower bounds for error bars
dat_cas$ub<-dat_cas$coef+1.96*dat_cas$se
dat_cas$lb<-dat_cas$coef-1.96*dat_cas$se

dat_cas$ub90<-dat_cas$coef+1.645*dat_cas$se
dat_cas$lb90<-dat_cas$coef-1.645*dat_cas$se

# Generate graph for Figure 4, casualty sensitivity
pdf(file="casualty_mechanism_new_garamond.pdf", width = 5.5, height = 3)
ggplot(dat_cas[c(1,4,7),], aes(names, coef, group = admin, color = admin))+
  geom_errorbar(ymin = dat_cas$lb[c(1,4,7)], ymax = dat_cas$ub[c(1,4,7)], 
                width = 0.2, 
                position=position_dodge(width=0.5))+
  geom_point(position=position_dodge(width=0.5), size = 0.8)+
  scale_y_continuous(name = "Regression Coefficient",limits = c(-0.5, 0.5),
                     breaks = c(-0.5,-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5))+
  geom_hline(yintercept = 0)+scale_x_discrete(name = "")+
  geom_vline(xintercept = 1.5, linetype = "longdash")+
  scale_x_discrete("", labels = c("IV: Black\nDV: Cas. Sens.\nSample: All",
                              "IV: Cas. Sens.\nDV: Support for War\nSample: Black",
                              "IV: Cas. Sens.\nDV: Support for War\nSample: White"))+
  scale_color_manual(name = "Presidency", values = c("grey55", "black"))+
  theme_bw()+theme_set(theme_bw(base_size = 10, base_family = "EB Garamond" ))
dev.off()

setEPS()                                             
postscript("casualty_mechanism_new_garamond.eps", width = 5.5, height = 3)
ggplot(dat_cas[c(1,4,7),], aes(names, coef, group = admin, color = admin))+
  geom_errorbar(ymin = dat_cas$lb[c(1,4,7)], ymax = dat_cas$ub[c(1,4,7)], 
                width = 0.2, 
                position=position_dodge(width=0.5))+
  geom_point(position=position_dodge(width=0.5), size = 0.8)+
  scale_y_continuous(name = "Regression Coefficient",limits = c(-0.5, 0.5),
                     breaks = c(-0.5,-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5))+
  geom_hline(yintercept = 0)+scale_x_discrete(name = "")+
  geom_vline(xintercept = 1.5, linetype = "longdash")+
  scale_x_discrete("", labels = c("IV: Black\nDV: Cas. Sens.\nSample: All",
                                  "IV: Cas. Sens.\nDV: Support for War\nSample: Black",
                                  "IV: Cas. Sens.\nDV: Support for War\nSample: White"))+
  scale_color_manual(name = "Presidency", values = c("grey55", "black"))+
  theme_bw()+theme_set(theme_bw(base_size = 10, base_family = "EB Garamond" ))
dev.off()

# Generate graph for Figure 6, military family
pdf(file="service_mechanism_new_garamond.pdf", width = 5.5, height = 3)
ggplot(dat_cas[-c(1,4,7),], aes(names, coef, 
                                group = admin, color = admin))+
  geom_errorbar(ymin = dat_cas$lb[-c(1,4,7)], ymax = dat_cas$ub[-c(1,4,7)], 
                width = 0.2, 
                position=position_dodge(width=0.5))+
  geom_point(position=position_dodge(width=0.5), size = 0.8)+
  scale_y_continuous(name = "Regression Coefficient",limits = c(-0.5, 0.5),
                     breaks = c(-0.5,-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5))+
  geom_hline(yintercept = 0)+scale_x_discrete(name = "")+
  geom_vline(xintercept = 1.5, linetype = "longdash")+
  scale_x_discrete("", labels = c("IV: Black\nDV: Mil. Family\nSample: All",
                                  "IV: Mil. Family\nDV: Support for War\nSample: Black",
                                  "IV: Mil. Family\nDV: Support for War\nSample: White"))+
  scale_color_manual(name = "Presidency", values = c("grey55", "black"))+
  theme_bw()+theme_set(theme_bw(base_size = 10, base_family = "EB Garamond" ))
dev.off()

setEPS()                                             
postscript("service_mechanism_new_garamond.eps", width = 5.5, height = 3)
ggplot(dat_cas[-c(1,4,7),], aes(names, coef, 
                                group = admin, color = admin))+
  geom_errorbar(ymin = dat_cas$lb[-c(1,4,7)], ymax = dat_cas$ub[-c(1,4,7)], 
                width = 0.2, 
                position=position_dodge(width=0.5))+
  geom_point(position=position_dodge(width=0.5), size = 0.8)+
  scale_y_continuous(name = "Regression Coefficient",limits = c(-0.5, 0.5),
                     breaks = c(-0.5,-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5))+
  geom_hline(yintercept = 0)+scale_x_discrete(name = "")+
  geom_vline(xintercept = 1.5, linetype = "longdash")+
  scale_x_discrete("", labels = c("IV: Black\nDV: Mil. Family\nSample: All",
                                  "IV: Mil. Family\nDV: Support for War\nSample: Black",
                                  "IV: Mil. Family\nDV: Support for War\nSample: White"))+
  scale_color_manual(name = "Presidency", values = c("grey55", "black"))+
  theme_bw()+theme_set(theme_bw(base_size = 10, base_family = "EB Garamond" ))
dev.off()

#####################################################################################
# Output for Figure 7 (Coefficient Plot, Spending)
#####################################################################################

# Coefficients from each regression
effects_spend<-c(coef(bush_iq_spend_reg_1)[2], # bush, domestic spending
               coef(obama_iq_spend_reg_1)[2], # obama, domestic spending
               coef(bush_iq_spend_reg_2)[2], # bush, domestic spending, black sample
               coef(obama_iq_spend_reg_2)[2], # obama, domestic spending, black sample
               coef(bush_iq_spend_reg_3)[2], # bush, domestic spending, white sample
               coef(obama_iq_spend_reg_3)[2]) # obama, domestic spending, white sample

# Standard errors from each regression
se_spend<-c(bush_iq_spend_se_1[2], # bush, domestic spending
          obama_iq_spend_se_1[2], # obama, domestic spending
          bush_iq_spend_se_2[2], # bush, domestic spending, black sample
          obama_iq_spend_se_2[2], # obama, domestic spending, black sample
          bush_iq_spend_se_3[2], # bush, domestic spending, white sample
          obama_iq_spend_se_3[2]) # obama, domestic spending, white sample

# Labels for each dataset used
poll_names_spend<-c("3A\nDom. Spending\n(Bush, Iraq)", "3B\nDom. Spending\n(Obama, Iraq)",
                     "3A\nDom. Spending\n(Bush, Iraq)", "3B\nDom. Spending\n(Obama, Iraq)",
                    "3A\nDom. Spending\n(Bush, Iraq)", "3B\nDom. Spending\n(Obama, Iraq)")

# Category of each coefficient
effect_type_spend<-c(rep("Black (All)", 2),
                   rep("Spending (Black)",2),rep("Spending (White)",2))

# Administration
admin_type_spend<-c("Bush", "Obama",
                  "Bush", "Obama",
                  "Bush", "Obama")

# Combine to one data frame
dat_spend<-data.frame(names = effect_type_spend,
                    poll = poll_names_spend,
                    admin = admin_type_spend,
                    coef = effects_spend,
                    se = se_spend)

# Generate upper and lower bounds for error bars
dat_spend$ub<-dat_spend$coef+1.96*dat_spend$se
dat_spend$lb<-dat_spend$coef-1.96*dat_spend$se

dat_spend$ub90<-dat_spend$coef+1.645*dat_spend$se
dat_spend$lb90<-dat_spend$coef-1.645*dat_spend$se

# Generate graph for Figure 7, Spending Preferences
pdf(file="spending_mechanism_new_garamond.pdf", width = 5.5, height = 3)
ggplot(dat_spend, aes(names, coef, 
                                group = admin, color = admin))+
  geom_errorbar(ymin = dat_spend$lb, ymax = dat_spend$ub, 
                width = 0.2, 
                position=position_dodge(width=0.5))+
  geom_point(position=position_dodge(width=0.5), size = .8)+
  scale_y_continuous(name = "Regression Coefficient",limits = c(-0.5, 0.5),
                     breaks = c(-0.5,-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5))+
  geom_hline(yintercept = 0)+scale_x_discrete(name = "")+
  geom_vline(xintercept = 1.5, linetype = "longdash")+
  scale_x_discrete("", labels = c("IV: Black\nDV: Spending\nSample: All",
                                  "IV: Spending\nDV: Support for War\nSample: Black",
                                  "IV: Spending\nDV: Support for War\nSample: White"))+
  scale_color_manual(name = "Presidency", values = c("grey55", "black"))+
  theme_bw()+theme_set(theme_bw(base_size = 10, base_family = "EB Garamond" ))
dev.off()

setEPS()                                             
postscript("spending_mechanism_new_garamond.eps", width = 5.5, height = 3)
ggplot(dat_spend, aes(names, coef, 
                      group = admin, color = admin))+
  geom_errorbar(ymin = dat_spend$lb, ymax = dat_spend$ub, 
                width = 0.2, 
                position=position_dodge(width=0.5))+
  geom_point(position=position_dodge(width=0.5), size = .8)+
  scale_y_continuous(name = "Regression Coefficient",limits = c(-0.5, 0.5),
                     breaks = c(-0.5,-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5))+
  geom_hline(yintercept = 0)+scale_x_discrete(name = "")+
  geom_vline(xintercept = 1.5, linetype = "longdash")+
  scale_x_discrete("", labels = c("IV: Black\nDV: Spending\nSample: All",
                                  "IV: Spending\nDV: Support for War\nSample: Black",
                                  "IV: Spending\nDV: Support for War\nSample: White"))+
  scale_color_manual(name = "Presidency", values = c("grey55", "black"))+
  theme_bw()+theme_set(theme_bw(base_size = 10, base_family = "EB Garamond" ))
dev.off()

#####################################################################################
# Output for Figure 9 (Coefficient Plot, Alienation)
#####################################################################################

# Coefficients from each regression
effects_alien<-c(coef(Bush1alien_IQ)[2], # bush, alienation
                 coef(Obama1alien_IQ)[2], # obama, alienation
                 coef(Bush2alien_IQ)[2], # bush, alienation, black sample
                 coef(Obama2alien_IQ)[2], # obama, alienation, black sample
                 coef(Bush3alien_IQ)[2], # bush, alienation, white sample
                 coef(Obama3alien_IQ)[2]) # obama, alienation, white sample

# Standard errors from each regression
se_alien<-c(Bush1alien_IQ_se[2], # bush, alienation
            Obama1alien_IQ_se[2], # obama, alienation
            Bush2alien_IQ_se[2], # bush, alienation, black sample
            Obama2alien_IQ_se[2], # obama, alienation, black sample
            Bush3alien_IQ_se[2], # bush, alienation, white sample
            Obama3alien_IQ_se[2]) # obama, alienation, white sample

# Labels for each dataset used
poll_names_alien<-c("4A\nAlienation\n(Bush, Iraq)", "4B\nAlienation\n(Obama, Iraq)",
                    "4A\nAlienation\n(Bush, Iraq)", "4B\nAlienation\n(Obama, Iraq)",
                    "4A\nAlienation\n(Bush, Iraq)", "4B\nAlienation\n(Obama, Iraq)")

# Category of each coefficient
effect_type_alien<-c(rep("Black (All)", 2),
                     rep("Spending (Black)",2),rep("Spending (White)",2))
# Administration
admin_type_alien<-c("Bush", "Obama",
                    "Bush", "Obama",
                    "Bush", "Obama")

# Combine to one data frame
dat_alien<-data.frame(names = effect_type_alien,
                      poll = poll_names_alien,
                      admin = admin_type_alien,
                      coef = effects_alien,
                      se = se_alien)

# Generate upper and lower bounds for error bars
dat_alien$ub<-dat_alien$coef+1.96*dat_alien$se
dat_alien$lb<-dat_alien$coef-1.96*dat_alien$se

dat_alien$ub90<-dat_alien$coef+1.645*dat_alien$se
dat_alien$lb90<-dat_alien$coef-1.645*dat_alien$se

# Generate Figure 9
pdf(file="alienation_mechanism_new_garamond.pdf", width = 5.5, height = 3)
ggplot(dat_alien, aes(names, coef, 
                      group = admin, color = admin))+
  geom_errorbar(ymin = dat_alien$lb, ymax = dat_alien$ub, 
                width = 0.2, 
                position=position_dodge(width=0.5))+
  geom_point(position=position_dodge(width=0.5), size = 0.8)+
  scale_y_continuous(name = "Regression Coefficient",limits = c(-0.5, 0.5),
                     breaks = c(-0.5,-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5))+
  geom_hline(yintercept = 0)+scale_x_discrete(name = "")+
  geom_vline(xintercept = 1.5, linetype = "longdash")+
  scale_x_discrete("", labels = c("IV: Black\nDV: Alien.\nSample: All",
                                  "IV: Alien.\nDV: Support for War\nSample: Black",
                                  "IV: Alien.\nDV: Support for War\nSample: White"))+
  scale_color_manual(name = "Presidency", values = c("grey55", "black"))+
  theme_bw()+theme_set(theme_bw(base_size = 10, base_family = "EB Garamond" ))
dev.off()

setEPS()                                             
postscript("alienation_mechanism_new_garamond.eps", width = 5.5, height = 3)
ggplot(dat_alien, aes(names, coef, 
                      group = admin, color = admin))+
  geom_errorbar(ymin = dat_alien$lb, ymax = dat_alien$ub, 
                width = 0.2, 
                position=position_dodge(width=0.5))+
  geom_point(position=position_dodge(width=0.5), size = 0.8)+
  scale_y_continuous(name = "Regression Coefficient",limits = c(-0.5, 0.5),
                     breaks = c(-0.5,-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5))+
  geom_hline(yintercept = 0)+scale_x_discrete(name = "")+
  geom_vline(xintercept = 1.5, linetype = "longdash")+
  scale_x_discrete("", labels = c("IV: Black\nDV: Alien.\nSample: All",
                                  "IV: Alien.\nDV: Support for War\nSample: Black",
                                  "IV: Alien.\nDV: Support for War\nSample: White"))+
  scale_color_manual(name = "Presidency", values = c("grey55", "black"))+
  theme_bw()+theme_set(theme_bw(base_size = 10, base_family = "EB Garamond" ))
dev.off()

# ####################################################
# # Output for Old Combined Coefficient Plot, "Black" (Various Hypotheses)
# ####################################################
# 
# # Coefficients from each mdoel for "black"
# coef_black_effects<-c(coef(bush_iq_cas_reg_1)[2], # bush, casualty sens
#               coef(bush_iq_serv_reg_1)[2], # bush, family member serving
#               coef(obama_iq_serv_reg_1)[2], # obama, family member serving
#               coef(bush_iq_spend_reg_1)[2],  # bush spending
#               coef(obama_iq_spend_reg_1)[2], # obama spending
#               coef(Bush1alien_IQ)[2], # bush alienation
#               coef(Obama1alien_IQ)[2]) # obama alienation
# 
# # SE from each model for "black"
# se_black_effects<-c(bush_iq_cas_se_1[2], # bush, casualty sens
#                     bush_iq_serv_se_1[2], # bush, family member serving
#                     obama_iq_serv_se_1[2], # obama, family member serving
#                     bush_iq_spend_se_1[2],  # bush spending
#                     obama_iq_spend_se_1[2], # obama spending
#                     Bush1alien_IQ_se[2], # bush alienation
#                     Obama1alien_IQ_se[2]) # obama alienation
# 
# # Labels for each dataset used
# poll_names_hyp<-c("1A\nCas. Sensitivity\n(Bush, Iraq)", 
#                   "2A\nMil. Family\n(Bush, Iraq)", "2B\nMil. Family\n(Obama, Iraq)", 
#                   "3A\nDom. Spending\n(Bush, Iraq)", "3B\nDom. Spending\n(Obama, Iraq)",
#                   "4A\nAlienation\n(Bush, Iraq)", "4B\nAlienation\n(Obama, Iraq)")
# 
# # Identifiers
# poll_type_hyp<-c("1. Casualty Sensitivity", 
#              "2. Military Family Member", "2. Military Family Member",
#              "3. Prefer Domestic Spending", "3. Prefer Domestic Spending",
#              "4. Political Alienation", "4. Political Alienation")
# 
# # Combine into data frame for visualization
# dat<-data.frame(names = poll_names_hyp,
#                 type = poll_type_hyp,
#                 coef = coef_black_effects,
#                 se = se_black_effects)
# 
# # Generate upper and lower bounds for error bars
# dat$ub<-dat$coef+1.96*dat$se
# dat$lb<-dat$coef-1.96*dat$se
# 
# dat$ub90<-dat$coef+1.645*dat$se
# dat$lb90<-dat$coef-1.645*dat$se
# 
# # Order the polls
# dat$names=factor(dat$names, levels = c("1A\nCas. Sensitivity\n(Bush, Iraq)", 
#                                        "2A\nMil. Family\n(Bush, Iraq)", "2B\nMil. Family\n(Obama, Iraq)", 
#                                        "3A\nDom. Spending\n(Bush, Iraq)", "3B\nDom. Spending\n(Obama, Iraq)",
#                                        "4A\nAlienation\n(Bush, Iraq)", "4B\nAlienation\n(Obama, Iraq)"))
# 
# # Generate graph for publication
# pdf(file="black_explanatory_vars.pdf", width = 10, height = 5)
# ggplot(dat, aes(names, coef, group = type, color = type))+
#   geom_errorbar(ymin = dat$lb, ymax = dat$ub, width = 0.2, 
#                 position=position_dodge(width=0.5))+
#   geom_point(position=position_dodge(width=0.5), size = 0.8)+
#   geom_point(aes(shape = type))+ 
#   scale_y_continuous(name = "Coefficient on 'Black'",limits = c(-0.25, 0.2),
#                      breaks = c(-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2))+
#   geom_hline(yintercept = 0)+scale_x_discrete(name = "")+
#   scale_shape_manual(name = "Outcome", values = c(1,2,8,4))+scale_color_manual(name = "Outcome", values = c("black", "grey55", "black", "grey55"))+theme_bw()
# dev.off()


####################################################
#
# 9) Repeats sections 2, 4-7 for Afghanistan War
#
####################################################

####################################################
# Regressions for Table 23, Figure 12 (Race Gap, Afghanistan War)
####################################################

# Models 1 & 2: Between-group & with controls, ABC polling dataset
bush_afgh_casualty_reg_1<-lm(support_war ~ black, 
                             data = poll_afgh)

bush_afgh_casualty_se_1<-sqrt(diag(vcovHC(bush_afgh_casualty_reg_1, type = "HC1")))

bush_afgh_casualty_reg_2<-lm(support_war ~ black + gender + 
                               party_id + education + age, 
                             data = poll_afgh)

bush_afgh_casualty_se_2<-sqrt(diag(vcovHC(bush_afgh_casualty_reg_2, type = "HC1")))

# Models 3 & 4: Between-group & with controls, Asstd. (Alienation) polling dataset, Bush
bush_afgh_reg_3<-lm(outcome ~ black+factor(poll), 
                    data = Bush_data_BW)
bush_afgh_se_3<-sqrt(diag(vcovHC(bush_afgh_reg_3, type = "HC1")))

bush_afgh_reg_4<-lm(outcome ~ black+gender+party_id+education+age+income+factor(poll), 
                    data = Bush_data_BW)
bush_afgh_se_4<-sqrt(diag(vcovHC(bush_afgh_reg_4, type = "HC1")))

# Models 5 & 6: Between-group & with controls, Asstd. (Alienation) polling dataset, Obama
obama_afgh_reg_5<-lm(outcome ~ black+factor(poll), 
                     data = Obama_data_BW)
obama_afgh_se_5<-sqrt(diag(vcovHC(obama_afgh_reg_5, type = "HC1")))

obama_afgh_reg_6<-lm(outcome ~ black+gender+party_id+education+age+income+factor(poll), 
                     data = Obama_data_BW)
obama_afgh_se_6<-sqrt(diag(vcovHC(obama_afgh_reg_6, type = "HC1")))


# Models 7 & 8: Between-group & with controls, CCES polling dataset, Obama
afgh_obama_cces_reg_9<-lm(afghanistan_support ~ black + factor(poll), data = cces_obama)

afgh_obama_cces_se_9<-sqrt(diag(vcovHC(afgh_obama_cces_reg_9,type = "HC1")))

afgh_obama_cces_reg_10<-lm(afghanistan_support ~ black + 
                             gender + party_id + education + 
                             age + income + factor(poll), data = cces_obama)

afgh_obama_cces_se_10<-sqrt(diag(vcovHC(afgh_obama_cces_reg_10,type = "HC1")))

####################################################
# Output for Table 23 (Regression Models 1-10)
####################################################
stargazer(bush_afgh_casualty_reg_1, bush_afgh_casualty_reg_2, 
          bush_afgh_reg_3, bush_afgh_reg_4, 
          obama_afgh_reg_5, obama_afgh_reg_6,
          afgh_obama_cces_reg_9, afgh_obama_cces_reg_10,
          se = list(bush_afgh_casualty_se_1, bush_afgh_casualty_se_2, 
                    bush_afgh_se_3, bush_afgh_se_4, 
                    obama_afgh_se_5, obama_afgh_se_6,
                    afgh_obama_cces_se_9, afgh_obama_cces_se_10),
          covariate.labels = c("Black", "Male",
                               "Independent", "Republican", "Any College",
                               "Some College", "College Grad", "Post-Grad",
                               "Age: 30-49", "Age: 50-64", "Age: 65+",
                               "Income: $30K-$50K",
                               "Income: $50K-$75K",
                               "Income: $75K and above",
                               rep(NA, 14), "Constant"),
          column.labels = c("Bush (Casualties)", "Bush (Alienation)",
                              "Obama (Alienation)", "Obama (CCES)"),
          column.separate = c(2,2,2,2),
          omit = c("poll"),
          omit.labels = c("Poll FE"),
          add.lines = list(c("Polls", 1,1,4,4,11,11,2,2)),
          omit.stat = c("adj.rsq","f", "ser"),
          out = "table_C1_afgh_full_new.tex")

# 14
####################################################
# Output for Figure 12 (Coefficient Plot, "Black," Models 1-10)
####################################################

# Coefficients from each mdoel for "black"
coef_black<-c(coef(bush_afgh_casualty_reg_1)[2],
              coef(bush_afgh_casualty_reg_2)[2],
              coef(bush_afgh_reg_3)[2],
              coef(bush_afgh_reg_4)[2],  
              coef(obama_afgh_reg_5)[2], 
              coef(obama_afgh_reg_6)[2], 
              coef(afgh_obama_cces_reg_9)[2],
              coef(afgh_obama_cces_reg_10)[2])

# SE from each model for "black"
se_black<-c(bush_afgh_casualty_se_1[2],
            bush_afgh_casualty_se_2[2], 
            bush_afgh_se_3[2],
            bush_afgh_se_4[2], 
            obama_afgh_se_5[2],
            obama_afgh_se_6[2], 
            afgh_obama_cces_se_9[2],
            afgh_obama_cces_se_10[2])

# Labels for each dataset used
poll_names<-c("ABC\n(Bush, Afgh.)", "ABC\n(Bush, Afgh.)",
              "Asstd.\n(Bush, Afgh.)", "Asstd.\n(Bush, Afgh.)", 
              "Asstd.\n(Obama, Afgh.)", "Asstd.\n(Obama, Afgh.)",
              "CCES\n(Obama, Afgh.)","CCES\n(Obama, Afgh.)")

# Identifiers
poll_type<-c("Bivariate", "Full Regression",
             "Bivariate", "Full Regression",
             "Bivariate", "Full Regression",
             "Bivariate", "Full Regression")

# Combine into data frame for visualization
dat<-data.frame(names = poll_names,
                type = poll_type,
                coef = coef_black,
                se = se_black)

# Generate upper and lower bounds for error bars
dat$ub<-dat$coef+1.96*dat$se
dat$lb<-dat$coef-1.96*dat$se

dat$ub90<-dat$coef+1.645*dat$se
dat$lb90<-dat$coef-1.645*dat$se

# Order the polls
dat$names=factor(dat$names, levels = c("ABC\n(Bush, Afgh.)", 
                                       "Asstd.\n(Bush, Afgh.)", 
                                       "Asstd.\n(Obama, Afgh.)",
                                       "CCES\n(Obama, Afgh.)"))

# Generate graph for publication
pdf(file="black_coef_afgh_new.pdf", width = 10, height = 5)
ggplot(dat, aes(names, coef, group = type, color = type))+
  geom_errorbar(ymin = dat$lb, ymax = dat$ub, width = 0.2, 
                position=position_dodge(width=0.5))+
  geom_point(position=position_dodge(width=0.5))+
  scale_y_continuous(name = "Coefficient on 'Black'",limits = c(-0.4, 0.2),
                     breaks = c(-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2))+
  geom_hline(yintercept = 0)+scale_x_discrete(name = "")+
  scale_color_manual(name = "Specification", values = c("grey55", "black"))+theme_bw()
dev.off()

###########################################
# Output for Table 24 - As above, Democrats only
###########################################

# Models 1 & 2: Between-group & with controls, ABC polling dataset
bush_afgh_casualty_dem_reg_1<-lm(support_war ~ black, 
                                 data = poll_afgh[which(poll_afgh$party_id=="democrat"),])

bush_afgh_casualty_dem_se_1<-sqrt(diag(vcovHC(bush_afgh_casualty_dem_reg_1, type = "HC1")))

bush_afgh_casualty_dem_reg_2<-lm(support_war ~ black + gender + 
                                   education + age, 
                                 data = poll_afgh[which(poll_afgh$party_id=="democrat"),])

bush_afgh_casualty_dem_se_2<-sqrt(diag(vcovHC(bush_afgh_casualty_dem_reg_2, type = "HC1")))

# Models 3 & 4: Between-group & with controls, Asstd. polling dataset, Bush
bush_afgh_dem_reg_3<-lm(outcome ~ black+factor(poll), 
                        data = Bush_data_BW[which(Bush_data_BW$party_id=="democrat"),])
bush_afgh_dem_se_3<-sqrt(diag(vcovHC(bush_afgh_dem_reg_3, type = "HC1")))

bush_afgh_dem_reg_4<-lm(outcome ~ black+gender+education+age+income+factor(poll), 
                        data = Bush_data_BW[which(Bush_data_BW$party_id=="democrat"),])
bush_afgh_dem_se_4<-sqrt(diag(vcovHC(bush_afgh_dem_reg_4, type = "HC1")))

# Models 5 & 6: Between-group & with controls, Asstd. polling dataset, Obama
obama_afgh_dem_reg_5<-lm(outcome ~ black+factor(poll), 
                         data = Obama_data_BW[which(Obama_data_BW$party_id=="democrat"),])
obama_afgh_dem_se_5<-sqrt(diag(vcovHC(obama_afgh_dem_reg_5, type = "HC1")))

obama_afgh_dem_reg_6<-lm(outcome ~ black+gender+education+age+income+factor(poll), 
                         data = Obama_data_BW[which(Obama_data_BW$party_id=="democrat"),])
obama_afgh_dem_se_6<-sqrt(diag(vcovHC(obama_afgh_dem_reg_6, type = "HC1")))

# Models 9 & 10: Between-group & with controls, CCES (Obama)
afgh_obama_cces_dem_reg_9<-lm(afghanistan_support ~ black + factor(poll), 
                              data = cces_obama[which(cces_obama$party_id=="democrat"),])

afgh_obama_cces_dem_se_9<-sqrt(diag(vcovHC(afgh_obama_cces_dem_reg_9,type = "HC1")))

afgh_obama_cces_dem_reg_10<-lm(afghanistan_support ~ black + 
                                 gender + education + 
                                 age + income + factor(poll), data = cces_obama[which(cces_obama$party_id=="democrat"),])

afgh_obama_cces_dem_se_10<-sqrt(diag(vcovHC(afgh_obama_cces_dem_reg_10,type = "HC1")))

##############################################################
# Output for Table 24 in the Appendix (Regression Models 1-10)
##############################################################
stargazer(bush_afgh_casualty_dem_reg_1, bush_afgh_casualty_dem_reg_2, 
          bush_afgh_dem_reg_3, bush_afgh_dem_reg_4, 
          obama_afgh_dem_reg_5, obama_afgh_dem_reg_6, 
          afgh_obama_cces_dem_reg_9, afgh_obama_cces_dem_reg_10,
          se = list(bush_afgh_casualty_dem_se_1, bush_afgh_casualty_dem_se_2, 
                    bush_afgh_dem_se_3, bush_afgh_dem_se_4, 
                    obama_afgh_dem_se_5, obama_afgh_dem_se_6,
                    afgh_obama_cces_dem_se_9, afgh_obama_cces_dem_se_10),
          covariate.labels = c("Black", "Male", "Any College",
                               "Some College", "College Grad", "Post-Grad",
                               "Age: 30-49", "Age: 50-64", "Age: 65+",
                               "Income: $30K-$50K",
                               "Income: $50K-$75K",
                               "Income: $75K and above",
                               rep(NA, 14), "Constant"),
          column.labels = c("Bush (Casualties)", "Bush (Alienation)",
                            "Obama (Alienation)", "Obama (CCES)"),
          column.separate = c(2,2,2,2),
          omit = c("poll"),
          omit.labels = c("Poll FE"),
          add.lines = list(c("Polls", 1,1,4,4,11,11,2,2)),
          omit.stat = c("adj.rsq","f", "ser"),
          out = "table_C2_afgh_dem_new.tex")


####################################################
# Generate Figure 13 in the Appendix
####################################################
dem_coef_black<-c(coef(bush_afgh_casualty_dem_reg_1)[2],
              coef(bush_afgh_casualty_dem_reg_2)[2],
              coef(bush_afgh_dem_reg_3)[2],
              coef(bush_afgh_dem_reg_4)[2],  
              coef(obama_afgh_dem_reg_5)[2], 
              coef(obama_afgh_dem_reg_6)[2], 
              coef(afgh_obama_cces_dem_reg_9)[2],
              coef(afgh_obama_cces_dem_reg_10)[2])

dem_se_black<-c(bush_afgh_casualty_dem_se_1[2],
                bush_afgh_casualty_dem_se_2[2], 
                bush_afgh_dem_se_3[2],
                bush_afgh_dem_se_4[2], 
                obama_afgh_dem_se_5[2],
                obama_afgh_dem_se_6[2], 
                afgh_obama_cces_dem_se_9[2],
                afgh_obama_cces_dem_se_10[2])

poll_names<-c("ABC\n(Bush, Afgh.)", "ABC\n(Bush, Afgh.)",
              "Asstd.\n(Bush, Afgh.)", "Asstd.\n(Bush, Afgh.)", 
              "Asstd.\n(Obama, Afgh.)", "Asstd.\n(Obama, Afgh.)",
              "CCES\n(Obama, Afgh.)","CCES\n(Obama, Afgh.)")

poll_type<-c("Bivariate", "Full Regression",
             "Bivariate", "Full Regression",
             "Bivariate", "Full Regression",
             "Bivariate", "Full Regression")

dat<-data.frame(names = poll_names,
                type = poll_type,
                coef = dem_coef_black,
                se = dem_se_black)

dat$ub<-dat$coef+1.96*dat$se
dat$lb<-dat$coef-1.96*dat$se

dat$ub90<-dat$coef+1.645*dat$se
dat$lb90<-dat$coef-1.645*dat$se

dat$names=factor(dat$names, levels = c("ABC\n(Bush, Afgh.)", 
                                       "Asstd.\n(Bush, Afgh.)", 
                                       "Asstd.\n(Obama, Afgh.)", 
                                       "CCES\n(Obama, Afgh.)"))

pdf(file="black_coef_afgh_dem_new.pdf", width = 10, height = 5)
ggplot(dat, aes(names, coef, group = type, color = type))+
  geom_errorbar(ymin = dat$lb, ymax = dat$ub, width = 0.2, 
                position=position_dodge(width=0.5))+
  geom_point(position=position_dodge(width=0.5))+
  scale_y_continuous(name = "Coefficient on 'Black'",limits = c(-0.4, 0.2),
                     breaks = c(-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2))+
  geom_hline(yintercept = 0)+scale_x_discrete(name = "")+
  scale_color_manual(name = "Specification", values = c("dodgerblue", "blue4"))+theme_bw()
dev.off()

###################################################################
#
# 10) Mechanisms - Afghanistan
#
###################################################################

###################################################################
#
# 10a) Costs of War - Casualties & Military Service, Afghanistan
#
###################################################################


###################################################################
# Regressions for Table 25, Models 1-3 (Casualties - Direct Effects)
###################################################################

# Model 1 - correlates of casualty sensitivity
bush_afgh_cas_reg_1<-lm(casualty_sens ~ black + gender + 
                          party_id + education + age, 
                      data = poll_afgh)

bush_afgh_cas_se_1<-sqrt(diag(vcovHC(bush_afgh_cas_reg_1, type = "HC1")))
coeftest(bush_afgh_cas_reg_1, vcovHC(bush_afgh_cas_reg_1, type = "HC1"))

# Model 2 - support for war, black sample only
bush_afgh_cas_reg_2<-lm(support_war ~ casualty_sens + gender + 
                        education + party_id + age, 
                      data = poll_afgh[which(poll_afgh$black==1),])

bush_afgh_cas_se_2<-sqrt(diag(vcovHC(bush_afgh_cas_reg_2, type = "HC1")))
coeftest(bush_afgh_cas_reg_2, vcovHC(bush_afgh_cas_reg_2, type = "HC1"))

# Model 3 - support for war, white sample only
bush_afgh_cas_reg_3<-lm(support_war ~ casualty_sens + gender + 
                          party_id + education + age, 
                      data = poll_afgh[which(poll_afgh$black==0),])

bush_afgh_cas_se_3<-sqrt(diag(vcovHC(bush_afgh_cas_reg_3, type = "HC1")))
coeftest(bush_afgh_cas_reg_3, vcovHC(bush_afgh_cas_reg_3, type = "HC1"))

# Mediation Analysis - Casualty Sensitivity
abc_med_afgh <-poll_afgh[which(!is.na(poll_afgh$support_war)),]

set.seed(02138)

med.fit.bush.cas.afgh <- lm(casualty_sens ~ black + gender +
                         education + party_id+ age, data=abc_med_afgh)
out.fit.bush.cas.afgh <- lm(support_war ~ black + casualty_sens + gender +
                         education + party_id+ age, data=abc_med_afgh)
med.out.bush.cas.afgh <- mediate(med.fit.bush.cas.afgh, out.fit.bush.cas.afgh,
                            treat = "black", mediator = "casualty_sens",
                            robustSE = TRUE, sims = 250)
summary(med.out.bush.cas.afgh)
plot(med.out.bush.cas.afgh, xlim = c(-0.35, 0.2))

###################################################################
# Regressions for Table 25, Models 4-6 (Military Service, Obama)
###################################################################

# Model 4 - Correlates of military service
obama_afgh_serv_reg_1<-lm(serving ~ black + 
                          gender +party_id+ education +
                          age +income + factor(poll), data = cces_obama[which(!is.na(cces_obama$afghanistan_support)),])

coeftest(obama_afgh_serv_reg_1, vcovHC(obama_afgh_serv_reg_1, 
                                     type = "HC1"))

obama_afgh_serv_se_1<-sqrt(diag(vcovHC(obama_afgh_serv_reg_1, 
                                     type = "HC1")))

# Model 5 - Correlates of support for war, black sample
obama_afgh_serv_reg_2<-lm(afghanistan_support ~ serving +
                          gender +party_id+ education +
                          age +income + factor(poll), 
                          data = cces_obama[which(cces_obama$black_aa==1),])

coeftest(obama_afgh_serv_reg_2, vcovHC(obama_afgh_serv_reg_2, 
                                     type = "HC1"))

obama_afgh_serv_se_2<-sqrt(diag(vcovHC(obama_afgh_serv_reg_2, 
                                     type = "HC1")))

# Model 6 - Correlates of support for war, white sample
obama_afgh_serv_reg_3<-lm(afghanistan_support ~ serving +
                          gender +party_id+ education +
                          age +income + factor(poll), 
                        data = cces_obama[which(cces_obama$black_aa==0),])

coeftest(obama_afgh_serv_reg_3, vcovHC(obama_afgh_serv_reg_3, 
                                     type = "HC1"))

obama_afgh_serv_se_3<-sqrt(diag(vcovHC(obama_afgh_serv_reg_3, 
                                     type = "HC1")))

# Mediation Analysis - Military Service, Obama
set.seed(02138)

med.fit.obama.serv.afgh <- lm(serving ~ black + 
                           gender + party_id + education +
                           age + income + factor(poll), data = cces_obama[which(!is.na(cces_obama$afghanistan_support)),])
out.fit.obama.serv.afgh <- lm(afghanistan_support ~ black + serving + 
                           gender + party_id + education +
                           age + income + factor(poll), data = cces_obama[which(!is.na(cces_obama$afghanistan_support)),])
med.out.obama.serv.afgh <- mediate(med.fit.obama.serv.afgh, out.fit.obama.serv.afgh, 
                              treat = "black", mediator = "serving",
                              robustSE = TRUE, sims = 250)
summary(med.out.obama.serv.afgh)
plot(med.out.obama.serv.afgh, xlim = c(-0.35, 0.2))

###################################################################
# Figure 14, Mediation Plot for Costs of War
###################################################################
pdf("casualty_mediation_plots_afgh_new.pdf", width = 4, height = 5)
par(mar=c(2, 4, 4, 2),
    mfcol = c(2, 1))
plot(med.out.bush.cas.afgh, xlim = c(-0.35, 0.2))
title(main = list("Casualty Sensitivity (Bush)", cex = 1.5,
                  col = "black", font = 1), line = 1, adj = 0)
plot(med.out.obama.serv.afgh, xlim = c(-0.35, 0.2))
title(main = list("Military Family (Obama)", cex = 1.5,
                  col = "black", font = 1), line = 1, adj = 0)
dev.off()

###################################################################
# Generates Table 25, Models 1-6 (Costs of War)
###################################################################
stargazer(bush_afgh_cas_reg_1,
          bush_afgh_cas_reg_2,
          bush_afgh_cas_reg_3,
          obama_afgh_serv_reg_1,
          obama_afgh_serv_reg_2,
          obama_afgh_serv_reg_3,
          se = list(bush_afgh_cas_se_1,
                    bush_afgh_cas_se_2,
                    bush_afgh_cas_se_3,
                    obama_afgh_serv_se_1,
                    obama_afgh_serv_se_2,
                    obama_afgh_serv_se_3),
          covariate.labels = c("Black", "Cas. Sens", "Mil. Family",
                               "Male",
                               "Independent", "Republican",
                               "Any College",
                               "Some College", "College Grad", "Post-Grad",
                               "Age: 30-49", "Age: 50-64", "Age: 65+",
                               "Income: $30K-$50K",
                               "Income: $50K-$75K",
                               "Income: $75K and above",
                               NA, "Constant"),
          column.labels = rep(c("Full", "Black", "White"), 2),
          omit = c("poll"),
          omit.labels = c("Poll FE"),
          add.lines = list(c("Polls", 1,1,1,2,2,2)),
          omit.stat = c("adj.rsq","f", "ser"),
          out = "table_C3_afgh_casualty_new.tex")

###################################################################
#
# 10b) Costs of War - Spending
#
###################################################################


###################################################################
# Regressions for Table 26, Models 1-3 (Spending, Obama)
###################################################################

# Model 1 - correlates of support for military spending, Obama
obama_afgh_spend_reg_1<-lm(dom_over_def ~ black + 
                           gender +party_id+ education +
                           age +income + factor(poll), data = cces_obama[which(!is.na(cces_obama$afghanistan_support)),])

coeftest(obama_afgh_spend_reg_1, vcovHC(obama_afgh_spend_reg_1, 
                                      type = "HC1"))

obama_afgh_spend_se_1<-sqrt(diag(vcovHC(obama_afgh_spend_reg_1, 
                                      type = "HC1")))

# Model 2 - correlates of support for war, black sample
obama_afgh_spend_reg_2<-lm(afghanistan_support ~ dom_over_def +
                           gender +party_id+ education +
                           age +income + factor(poll), data = cces_obama[which(cces_obama$black_aa==1),])

coeftest(obama_afgh_spend_reg_2, vcovHC(obama_afgh_spend_reg_2, 
                                      type = "HC1"))

obama_afgh_spend_se_2<-sqrt(diag(vcovHC(obama_afgh_spend_reg_2, 
                                      type = "HC1")))

# Model 3 - correlates of support for war, white sample
obama_afgh_spend_reg_3<-lm(afghanistan_support ~ dom_over_def +
                           gender +party_id+ education +
                           age +income + factor(poll), 
                         data = cces_obama[which(cces_obama$black_aa==0),])

coeftest(obama_afgh_spend_reg_3, vcovHC(obama_afgh_spend_reg_3, 
                                      type = "HC1"))

obama_afgh_spend_se_3<-sqrt(diag(vcovHC(obama_afgh_spend_reg_3, 
                                      type = "HC1")))

# Mediation Analysis - Spending, Obama
set.seed(02138)

med.fit.obama.spend.afgh <- lm(dom_over_def ~ black + 
                            gender + party_id + education +
                            age + income + factor(poll), data = cces_obama[which(!is.na(cces_obama$afghanistan_support)),])
out.fit.obama.spend.afgh <- lm(afghanistan_support ~ black + dom_over_def + 
                            gender + party_id + education +
                            age + income + factor(poll), data = cces_obama[which(!is.na(cces_obama$afghanistan_support)),])
med.out.obama.spend.afgh <- mediate(med.fit.obama.spend.afgh, out.fit.obama.spend.afgh, treat = "black", mediator = "dom_over_def",
                               robustSE = TRUE, sims = 250)
summary(med.out.obama.spend.afgh)
plot(med.out.obama.spend.afgh, xlim = c(-0.25, 0.2))

###################################################################
# Figure 15, Mediation Plot for Spending & War
###################################################################
pdf("spending_mediation_plots_afgh_new.pdf", width = 4, height = 2)
par(mar=c(2, 4, 4, 2),
    mfcol = c(1, 1))
plot(med.out.obama.spend.afgh, xlim = c(-0.2, 0.2), cex.axis = .8)
title(main = list("Spending (Obama)", cex = 1,
                  col = "black", font = 1), line = 1, adj = 0)

dev.off()

###################################################################
# Generates Table 26, Models 1-3 (Spending & War)
###################################################################
stargazer(obama_afgh_spend_reg_1,
          obama_afgh_spend_reg_2,
          obama_afgh_spend_reg_3,
          se = list(obama_afgh_spend_se_1,
                    obama_afgh_spend_se_2,
                    obama_afgh_spend_se_3),
          covariate.labels = c("Black", "Dom. Spending",
                               "Male",
                               "Independent", "Republican",
                               "Some College", "College Grad", "Post-Grad",
                               "Age: 30-49", "Age: 50-64", "Age: 65+",
                               "Income: $30K-$50K",
                               "Income: $50K-$75K",
                               "Income: $75K and above",
                               NA, "Constant"),
          column.labels = rep(c("Full", "Black", "White"), 1),
          omit = c("poll"),
          omit.labels = c("Poll FE"),
          add.lines = list(c("Polls", 2,2,2)),
          omit.stat = c("adj.rsq","f", "ser"),
          out = "table_C4_afgh_spending_new.tex")

###################################################################
#
# 10c) Prerequisites for War - Alienation
#
###################################################################

###################################################################
# Regressions for Table 27, Models 1-3 (Alienation, Bush)
###################################################################

# Model 1 - correlates of alienation, Bush
Bush1alien_afgh <- lm(alienation_binary ~ 
                      black+gender+party_id+education+age+income+factor(poll), data = Bush_data_BW)
Bush1alien_afgh_se<-sqrt(diag(vcovHC(Bush1alien_afgh, type = "HC1")))

# Model 2 - correlates of support for war, black sample
Bush2alien_afgh <- lm(outcome ~ alienation_binary +
                      gender+party_id+education+age+income+factor(poll), data = Bush_data_BW[which(Bush_data_BW$black==1),])
Bush2alien_afgh_se<-sqrt(diag(vcovHC(Bush2alien_afgh, type = "HC1")))

# Model 3 - correlates of support for war, white sample
Bush3alien_afgh <- lm(outcome ~ alienation_binary +
                      gender+party_id+education+age+income+factor(poll), data = Bush_data_BW[which(Bush_data_BW$black==0),])
Bush3alien_afgh_se<-sqrt(diag(vcovHC(Bush3alien_afgh, type = "HC1")))

# Mediation - Bush & Alienation
set.seed(02138)

med.fit.bush.alien.afgh <- lm(alienation_binary ~ black + 
                           gender+party_id+education+age+income+factor(poll), data = Bush_data_BW[!is.na(Bush_data_BW$outcome),])
out.fit.bush.alien.afgh <- lm(outcome ~ black + alienation_binary + 
                           gender+party_id+education+age+income+factor(poll), data = Bush_data_BW[!is.na(Bush_data_BW$outcome),])
med.out.bush.alien.afgh <- mediate(med.fit.bush.alien.afgh, out.fit.bush.alien.afgh, 
                              treat = "black", mediator = "alienation_binary",
                              robustSE = TRUE, sims = 250)
summary(med.out.bush.alien.afgh)
plot(med.out.bush.alien.afgh, xlim = c(-0.25, 0.2))


###################################################################
# Regressions for Table 27, Models 4-6 (Alienation, Obama)
###################################################################

# Model 4 - correlates of alienation, Obama
Obama1alien_afgh <- lm(alienation_binary ~ 
                       black+gender+party_id+education+age+income+factor(poll), data = Obama_data_BW)
Obama1alien_afgh_se<-sqrt(diag(vcovHC(Obama1alien_afgh, type = "HC1")))

# Model 5 - correlates of support for war, black sample
Obama2alien_afgh <- lm(outcome ~ alienation_binary +
                       gender+party_id+education+age+income+factor(poll), data = Obama_data_BW[which(Obama_data_BW$black==1),])
Obama2alien_afgh_se<-sqrt(diag(vcovHC(Obama2alien_afgh, type = "HC1")))

# Model 6 - correlates of support for war, white sample
Obama3alien_afgh <- lm(outcome ~ alienation_binary +
                       gender+party_id+education+age+income+factor(poll), data = Obama_data_BW[which(Obama_data_BW$black==0),])
Obama3alien_afgh_se<-sqrt(diag(vcovHC(Obama3alien_afgh, type = "HC1")))

summary(Obama3alien_afgh)

# Mediation - Obama & Alienation
set.seed(02138)

med.fit.obama.alien.afgh <- lm(alienation_binary ~ black + 
                            gender+party_id+education+age+income+factor(poll), 
                          data = Obama_data_BW[which(!is.na(Obama_data_BW$outcome)),])
out.fit.obama.alien.afgh <- lm(outcome ~ black + alienation_binary + 
                            gender+party_id+education+age+income+factor(poll), 
                          data = Obama_data_BW[which(!is.na(Obama_data_BW$outcome)),])
med.out.obama.alien.afgh <- mediate(med.fit.obama.alien.afgh, out.fit.obama.alien.afgh, 
                               treat = "black", mediator = "alienation_binary",
                               robustSE = TRUE, sims = 250)
summary(med.out.obama.alien.afgh)
plot(med.out.obama.alien.afgh, xlim = c(-0.25, 0.2))

###################################################################
# Generates Table 27, Models 1-6 (Alienation & War)
###################################################################
stargazer(Bush1alien_afgh,Bush2alien_afgh,Bush3alien_afgh,
          Obama1alien_afgh,Obama2alien_afgh,Obama3alien_afgh,
          se = list(Bush1alien_afgh_se,Bush2alien_afgh_se,Bush3alien_afgh_se,
                    Obama1alien_afgh_se,Obama2alien_afgh_se,Obama3alien_afgh_se),
          covariate.labels = c("Black", "Alienation",
                               "Male",
                               "Independent", "Republican",
                               "Some College", "College Grad", "Post-Grad",
                               "Age: 30-49", "Age: 50-64", "Age: 65+",
                               "Income: $30K-$50K",
                               "Income: $50K-$75K",
                               "Income: $75K and above",
                               rep(NA,13), "Constant"),
          column.labels = rep(c("Full", "Black", "White"), 2),
          omit = c("poll"),
          omit.labels = c("Poll FE"),
          add.lines = list(c("Polls", 4,4,4,11,11,11)),
          omit.stat = c("adj.rsq","f", "ser"),
          out = "table_C5_afgh_alienation_new.tex")

###################################################################
# Figure 16, Mediation Plot for Alienation & War
###################################################################
pdf("alien_mediation_plots_afgh_new.pdf", width = 4, height = 4)
par(mar=c(2, 4, 4, 2),
    mfcol = c(2, 1))
plot(med.out.bush.alien.afgh, xlim = c(-0.2, 0.2), cex.axis = .8)
title(main = list("Alienation (Bush)", cex = 1,
                  col = "black", font = 1), line = 1, adj = 0)
plot(med.out.obama.alien.afgh, xlim = c(-0.2, 0.2), cex.axis = .8)
title(main = list("Alienation (Obama)", cex = 1,
                  col = "black", font = 1), line = 1, adj = 0)

dev.off()

####################################################################
#
# 11) Robustness check: Missing Income
#
####################################################################

###################################################################
# 11a) Regressions for Table A19, Models 1-4 & Figure A4 (Casualties - Direct Effects)
###################################################################

# Model 0 - correlates of casualty sensitivity
bush_iq_cas_noincome_reg_0<-lm(support_war ~ black + gender + 
                                 party_id + education + age + 
                                 factor(poll), 
                               data = abc_iraq_combine)

bush_iq_cas_noincome_se_0<-sqrt(diag(vcovHC(bush_iq_cas_noincome_reg_0, type = "HC1")))
coeftest(bush_iq_cas_noincome_reg_0, vcovHC(bush_iq_cas_noincome_reg_0, type = "HC1"))

# Model 1 - correlates of casualty sensitivity
bush_iq_cas_noincome_reg_1<-lm(casual_sens ~ black + gender + 
                          party_id + education + age + 
                        factor(poll), 
                      data = abc_iraq_combine)

bush_iq_cas_noincome_se_1<-sqrt(diag(vcovHC(bush_iq_cas_noincome_reg_1, type = "HC1")))
coeftest(bush_iq_cas_noincome_reg_1, vcovHC(bush_iq_cas_noincome_reg_1, type = "HC1"))

# Model 2 - support for war, black sample only
bush_iq_cas_noincome_reg_2<-lm(support_war ~ casual_sens + gender + 
                        party_id + education + age + 
                        factor(poll), 
                      data = abc_iraq_combine[which(abc_iraq_combine$black==1),])

bush_iq_cas_noincome_se_2<-sqrt(diag(vcovHC(bush_iq_cas_noincome_reg_2, type = "HC1")))
coeftest(bush_iq_cas_noincome_reg_2, vcovHC(bush_iq_cas_noincome_reg_2, type = "HC1"))

# Model 3 - support for war, white sample only
bush_iq_cas_noincome_reg_3<-lm(support_war ~ casual_sens + gender + 
                        party_id + education + age + 
                        factor(poll), 
                      data = abc_iraq_combine[which(abc_iraq_combine$black==0),])

bush_iq_cas_noincome_se_3<-sqrt(diag(vcovHC(bush_iq_cas_noincome_reg_3, type = "HC1")))
coeftest(bush_iq_cas_noincome_reg_3, vcovHC(bush_iq_cas_noincome_reg_3, type = "HC1"))

# Mediation Analysis - Casualty Sensitivity
abc_med_noinc <- abc_iraq_combine[-which(is.na(abc_iraq_combine$support_war)),]

set.seed(02138)

med.fit.bush.cas.inc <- lm(casual_sens ~ black + gender +
                         education + party_id+ age + 
                         factor(poll), data=abc_med_noinc)
out.fit.bush.cas.inc <- lm(support_war ~ black + casual_sens + gender +
                         education + party_id+ age + 
                         factor(poll), data=abc_med_noinc)
med.out.bush.cas.inc <- mediate(med.fit.bush.cas.inc, out.fit.bush.cas.inc, treat = "black", mediator = "casual_sens",
                            robustSE = TRUE, sims = 250)
summary(med.out.bush.cas.inc)
plot(med.out.bush.cas.inc, xlim = c(-0.4, 0.1))


stargazer(bush_iq_cas_noincome_reg_0,bush_iq_cas_noincome_reg_1,
          bush_iq_cas_noincome_reg_2,bush_iq_cas_noincome_reg_3,
          se=list(bush_iq_cas_noincome_se_0,bush_iq_cas_noincome_se_1,
                  bush_iq_cas_noincome_se_2,bush_iq_cas_noincome_se_3),
          covariate.labels = c("Black", "Cas. Sensitivity",
                               "Male",
                               "Independent", "Republican",
                               "Some College", "College Grad", "Post-Grad",
                               "Age: 30-49", "Age: 50-64", "Age: 65+",
                               rep(NA,17), "Constant"),
          column.labels = c("Full", "Full", "Black", "White"),
          omit = c("poll"),
          omit.labels = c("Poll FE"),
          add.lines = list(c("Polls", 19,19,19,19)),
          omit.stat = c("adj.rsq","f", "ser"),
          out = "table_A19_cas_no_income_new.tex")

pdf("casualty_mediation_no_income_new.pdf", width = 4, height = 2)
par(mar=c(2, 4, 4, 2),
    mfcol = c(1, 1))
plot(med.out.bush.cas.inc, xlim = c(-0.25, 0.2), cex.axis = .8)
title(main = list("Casualties (Bush)", cex = 1,
                  col = "black", font = 1), line = 1, adj = 0)
dev.off()

###################################################################
# 11b) Regressions for Table A20, Models 1-4 & Figure A5 (Alienation, Bush)
###################################################################

# Model 0 - correlates of alienation, Bush
Bush0alien_noincome_IQ <- lm(outcome ~ black+gender+party_id+education+age+factor(poll), data = Bush_data_IQ_BW)
Bush0alien_noincome_IQ_se<-sqrt(diag(vcovHC(Bush0alien_noincome_IQ, type = "HC1")))

# Model 1 - correlates of alienation, Bush
Bush1alien_noincome_IQ <- lm(alienation_binary ~ 
                      black+gender+party_id+education+age+factor(poll), data = Bush_data_IQ_BW)
Bush1alien_noincome_IQ_se<-sqrt(diag(vcovHC(Bush1alien_noincome_IQ, type = "HC1")))

# Model 2 - correlates of support for war, black sample
Bush2alien_noincome_IQ <- lm(outcome ~ alienation_binary +
                      gender+party_id+education+age+factor(poll), data = Bush_data_IQ_BW[which(Bush_data_IQ_BW$black==1),])
Bush2alien_noincome_IQ_se<-sqrt(diag(vcovHC(Bush2alien_noincome_IQ, type = "HC1")))

# Model 3 - correlates of support for war, white sample
Bush3alien_noincome_IQ <- lm(outcome ~ alienation_binary +
                      gender+party_id+education+age+factor(poll), data = Bush_data_IQ_BW[which(Bush_data_IQ_BW$black==0),])
Bush3alien_noincome_IQ_se<-sqrt(diag(vcovHC(Bush3alien_noincome_IQ, type = "HC1")))

# Mediation - Bush & Alienation
set.seed(02138)

med.fit.bush.alien.inc <- lm(alienation_binary ~ black + 
                           gender+party_id+education+age+factor(poll), data = Bush_data_IQ_BW[!is.na(Bush_data_IQ_BW$outcome),])
out.fit.bush.alien.inc <- lm(outcome ~ black + alienation_binary + 
                           gender+party_id+education+age+factor(poll), data = Bush_data_IQ_BW[!is.na(Bush_data_IQ_BW$outcome),])
med.out.bush.alien.inc <- mediate(med.fit.bush.alien.inc, out.fit.bush.alien.inc, 
                              treat = "black", mediator = "alienation_binary",
                              robustSE = TRUE, sims = 250)
summary(med.out.bush.alien.inc)
plot(med.out.bush.alien.inc, xlim = c(-0.25, 0.2))

stargazer(Bush0alien_noincome_IQ, Bush1alien_noincome_IQ,
          Bush2alien_noincome_IQ, Bush3alien_noincome_IQ,
          se = list(Bush0alien_noincome_IQ_se,
                    Bush1alien_noincome_IQ_se,
                    Bush2alien_noincome_IQ_se,
                    Bush3alien_noincome_IQ_se),
          covariate.labels = c("Black", "Alienation",
                               "Male",
                               "Independent", "Republican",
                               "Some College", "College Grad", "Post-Grad",
                               "Age: 30-49", "Age: 50-64", "Age: 65+",
                               rep(NA,38), "Constant"),
          column.labels = c("Full", "Full", "Black", "White"),
          omit = c("poll"),
          omit.labels = c("Poll FE"),
          add.lines = list(c("Polls", 39,39,39,39)),
          omit.stat = c("adj.rsq","f", "ser"),
          out = "table_A20_alien_no_income_new.tex")

pdf("alienation_mediation_no_income_new.pdf", width = 4, height = 2)
par(mar=c(2, 4, 4, 2),
    mfcol = c(1, 1))
plot(med.out.bush.alien.inc, xlim = c(-0.2, 0.2), cex.axis = .8)
title(main = list("Alienation (Bush)", cex = 1,
                  col = "black", font = 1), line = 1, adj = 0)
dev.off()

####################################################################
#
# 12) Robustness check: Alternative Casualties Measure
#
####################################################################

###################################################################
# 12a) Regressions for Table A21, Models 1-3 & Figure A6 (Casualties - Alternative Measure)
###################################################################

# Model 1 - correlates of casualty sensitivity
bush_iq_cas_reg_1<-lm(casualty_unaccept ~ black + gender + 
                        party_id + education + age + income +
                        factor(poll), 
                      data = abc_iraq_combine)

bush_iq_cas_se_1<-sqrt(diag(vcovHC(bush_iq_cas_reg_1, type = "HC1")))
coeftest(bush_iq_cas_reg_1, vcovHC(bush_iq_cas_reg_1, type = "HC1"))

# Model 2 - support for war, black sample only
bush_iq_cas_reg_2<-lm(support_war ~ casualty_unaccept + gender + 
                        party_id + education + age + income +
                        factor(poll), 
                      data = abc_iraq_combine[which(abc_iraq_combine$black==1),])

bush_iq_cas_se_2<-sqrt(diag(vcovHC(bush_iq_cas_reg_2, type = "HC1")))
coeftest(bush_iq_cas_reg_2, vcovHC(bush_iq_cas_reg_2, type = "HC1"))

# Model 3 - support for war, white sample only
bush_iq_cas_reg_3<-lm(support_war ~ casualty_unaccept + gender + 
                        party_id + education + age + income +
                        factor(poll), 
                      data = abc_iraq_combine[which(abc_iraq_combine$black==0),])

bush_iq_cas_se_3<-sqrt(diag(vcovHC(bush_iq_cas_reg_3, type = "HC1")))
coeftest(bush_iq_cas_reg_3, vcovHC(bush_iq_cas_reg_3, type = "HC1"))

# Mediation Analysis - Casualty Sensitivity
abc_med_altcas <- abc_iraq_combine[-which(is.na(abc_iraq_combine$support_war)),]

set.seed(02138)

med.fit.bush.cas.alt <- lm(casualty_unaccept ~ black + gender +
                         education + party_id+ age + income +
                         factor(poll), data=abc_med_altcas)
out.fit.bush.cas.alt <- lm(support_war ~ black + casualty_unaccept + gender +
                         education + party_id+ age + income +
                         factor(poll), data=abc_med_altcas)
med.out.bush.cas.alt <- mediate(med.fit.bush.cas.alt, out.fit.bush.cas.alt, treat = "black", mediator = "casualty_unaccept",
                            robustSE = TRUE, sims = 250)
summary(med.out.bush.cas.alt)
plot(med.out.bush.cas.alt, xlim = c(-0.25, 0.2))

############################
stargazer(bush_iq_cas_reg_1,bush_iq_cas_reg_2,bush_iq_cas_reg_3,
          se = list(bush_iq_cas_se_1,bush_iq_cas_se_2,bush_iq_cas_se_3),
          covariate.labels = c("Black", "Casualty Sens.", "Male",
                               "Independent", "Republican",
                               "Some College", "College Grad", "Post-Grad",
                               "Age: 30-49", "Age: 50-64", "Age: 65+",
                               "Income: $30K-$50K",
                               "Income: $50K-$75K",
                               "Income: $75K and above",
                               rep(NA, 6), "Constant"),
          omit = c("poll"),
          column.labels = c("Full", "Black", "White"),
          omit.labels = c("Poll FE"),
          add.lines = list(c("Polls", 7,7,7)),
          omit.stat = c("adj.rsq","f", "ser"),
          out = "table_A21_cas_regression_new.tex")

pdf("casualty_mediation_alt_new.pdf", width = 4, height = 2)
par(mar=c(2, 4, 4, 2),
    mfcol = c(1, 1))
plot(med.out.bush.cas.alt, xlim = c(-0.25, 0.2), cex.axis = .8)
title(main = list("Casualties Alt. (Bush)", cex = 1,
                  col = "black", font = 1), line = 1, adj = 0)
dev.off()

###################################################################
#
# 13) Risks of War - Casualties & Military Service (Bivariate)
#
###################################################################

####################################################################################
# 13a) Regressions for Figure A1 & Table A16, Models 1-3 (Casualties - Direct Effects)
####################################################################################

# Model 1 - correlates of casualty sensitivity
bush_iq_cas_bin_reg_1<-lm(casual_sens ~ black+
                            factor(poll), 
                          data = abc_iraq_combine)

bush_iq_cas_bin_se_1<-sqrt(diag(vcovHC(bush_iq_cas_bin_reg_1, type = "HC1")))
coeftest(bush_iq_cas_bin_reg_1, vcovHC(bush_iq_cas_bin_reg_1, type = "HC1"))

# Model 2 - support for war, black sample only
bush_iq_cas_bin_reg_2<-lm(support_war ~ casual_sens+
                            factor(poll), 
                          data = abc_iraq_combine[which(abc_iraq_combine$black==1),])

bush_iq_cas_bin_se_2<-sqrt(diag(vcovHC(bush_iq_cas_bin_reg_2, type = "HC1")))
coeftest(bush_iq_cas_bin_reg_2, vcovHC(bush_iq_cas_bin_reg_2, type = "HC1"))

# values
coef(bush_iq_cas_bin_reg_1)[2]+1.96*bush_iq_cas_bin_se_1[2]
coef(bush_iq_cas_bin_reg_1)[2]-1.96*bush_iq_cas_bin_se_1[2]

# Model 3 - support for war, white sample only
bush_iq_cas_bin_reg_3<-lm(support_war ~ casual_sens+
                            factor(poll), 
                          data = abc_iraq_combine[which(abc_iraq_combine$black==0),])

bush_iq_cas_bin_se_3<-sqrt(diag(vcovHC(bush_iq_cas_bin_reg_3, type = "HC1")))
coeftest(bush_iq_cas_bin_reg_3, vcovHC(bush_iq_cas_bin_reg_3, type = "HC1"))

# Mediation Analysis - Casualty Sensitivity
abc_med_binary <- abc_iraq_combine[-which(is.na(abc_iraq_combine$support_war)),]

set.seed(02138)

med.fit.bush.cas.bin <- lm(casual_sens ~ black+
                             factor(poll), data=abc_med_binary)
out.fit.bush.cas.bin <- lm(support_war ~ black + casual_sens+
                             factor(poll), data=abc_med_binary)
med.out.bush.cas.bin <- mediate(med.fit.bush.cas.bin, out.fit.bush.cas.bin, treat = "black", mediator = "casual_sens",
                                robustSE = TRUE, sims = 250)
summary(med.out.bush.cas.bin)
plot(med.out.bush.cas.bin, xlim = c(-0.4, 0.1))

###############################################################################
# 13b) Regressions for Figure 6 & Table A16, Models 12-6 (Military Service, Bush, Bivariate)
###############################################################################

# Model 1 - Correlates of military service
bush_iq_serv_bin_reg_1<-lm(serving ~ black+
                             factor(poll), data = cces_bush)

coeftest(bush_iq_serv_bin_reg_1, vcovHC(bush_iq_serv_bin_reg_1, 
                                        type = "HC1"))

bush_iq_serv_bin_se_1<-sqrt(diag(vcovHC(bush_iq_serv_bin_reg_1, 
                                        type = "HC1")))

# Model 5 - Correlates of support for war, black sample
bush_iq_serv_bin_reg_2<-lm(iraq_support ~ serving+
                             factor(poll), data = cces_bush[which(cces_bush$black_aa==1),])

coeftest(bush_iq_serv_bin_reg_2, vcovHC(bush_iq_serv_bin_reg_2, 
                                        type = "HC1"))

bush_iq_serv_bin_se_2<-sqrt(diag(vcovHC(bush_iq_serv_bin_reg_2, 
                                        type = "HC1")))
# Model 6 - Correlates of support for war, white sample
bush_iq_serv_bin_reg_3<-lm(iraq_support ~ serving+
                             factor(poll), data = cces_bush[which(cces_bush$black_aa==0),])

coeftest(bush_iq_serv_bin_reg_3, vcovHC(bush_iq_serv_bin_reg_3, 
                                        type = "HC1"))

bush_iq_serv_bin_se_3<-sqrt(diag(vcovHC(bush_iq_serv_bin_reg_3, 
                                        type = "HC1")))


# Mediation Analysis - Military Service, Bush
set.seed(02138)

med.fit.bush.serv.bin <- lm(serving ~ black+factor(poll), data = cces_bush[-which(is.na(cces_bush$iraq_support)),])
out.fit.bush.serv.bin <- lm(iraq_support ~ black + serving+
                              factor(poll), data = cces_bush[-which(is.na(cces_bush$iraq_support)),])
med.out.bush.serv.bin <- mediate(med.fit.bush.serv.bin, out.fit.bush.serv.bin, treat = "black", mediator = "serving",
                                 robustSE = TRUE, sims = 250)
summary(med.out.bush.serv.bin)
plot(med.out.bush.serv.bin, xlim = c(-0.4, 0.1))

###################################################################
# 13c) Regressions for Figure 6 & Table A4, Models 7-9 (Military Service, Obama, Bivariate)
###################################################################

# Model 7 - Correlates of military service
obama_iq_serv_bin_reg_1<-lm(serving ~ black+
                              factor(poll), data = cces_obama)

coeftest(obama_iq_serv_bin_reg_1, vcovHC(obama_iq_serv_bin_reg_1, 
                                         type = "HC1"))

obama_iq_serv_bin_se_1<-sqrt(diag(vcovHC(obama_iq_serv_bin_reg_1, 
                                         type = "HC1")))

# Model 8 - Correlates of support for war, black sample
obama_iq_serv_bin_reg_2<-lm(iraq_support ~ serving+
                              factor(poll), 
                            data = cces_obama[which(cces_obama$black_aa==1),])

coeftest(obama_iq_serv_bin_reg_2, vcovHC(obama_iq_serv_bin_reg_2, 
                                         type = "HC1"))

obama_iq_serv_bin_se_2<-sqrt(diag(vcovHC(obama_iq_serv_bin_reg_2, 
                                         type = "HC1")))

# Model 9 - Correlates of support for war, white sample
obama_iq_serv_bin_reg_3<-lm(iraq_support ~ serving+
                              factor(poll), 
                            data = cces_obama[which(cces_obama$black_aa==0),])

coeftest(obama_iq_serv_bin_reg_3, vcovHC(obama_iq_serv_bin_reg_3, 
                                         type = "HC1"))

obama_iq_serv_bin_se_3<-sqrt(diag(vcovHC(obama_iq_serv_bin_reg_3, 
                                         type = "HC1")))

# Mediation Analysis - Military Service, Obama
set.seed(02138)

med.fit.obama.serv.bin <- lm(serving ~ black + factor(poll), data = cces_obama[-which(is.na(cces_obama$iraq_support)),])
out.fit.obama.serv.bin <- lm(iraq_support ~ black + serving + factor(poll), data = cces_obama[-which(is.na(cces_obama$iraq_support)),])
med.out.obama.serv.bin <- mediate(med.fit.obama.serv.bin, out.fit.obama.serv.bin, 
                                  treat = "black", mediator = "serving",
                                  robustSE = TRUE, sims = 250)
summary(med.out.obama.serv.bin)
plot(med.out.obama.serv.bin, xlim = c(-0.4, 0.1))

###################################################################
# 13d) Figure A1, Mediation Plot for Costs of War, Bivariate
###################################################################
pdf("casualty_bin_mediation_plots_new.pdf", width = 4, height = 5)
par(mar=c(2, 4, 4, 2),
    mfcol = c(3,1))
plot(med.out.bush.cas.bin, xlim = c(-0.4, 0.1))
title(main = list("Casualty Sensitivity (Bush)", cex = 1.5,
                  col = "black", font = 1), line = 1, adj = 0)
plot(med.out.bush.serv.bin, xlim = c(-0.4, 0.1))
title(main = list("Military Family (Bush)", cex = 1.5,
                  col = "black", font = 1), line = 1, adj = 0)
plot(med.out.obama.serv.bin, xlim = c(-0.4, 0.1))
title(main = list("Military Family (Obama)", cex = 1.5,
                  col = "black", font = 1), line = 1, adj = 0)
dev.off()

###################################################################
# 13e) Generates Table A16, Models 1-9 (Costs of War)
###################################################################
stargazer(bush_iq_cas_bin_reg_1,
          bush_iq_cas_bin_reg_2,
          bush_iq_cas_bin_reg_3,
          bush_iq_serv_bin_reg_1,
          bush_iq_serv_bin_reg_2,
          bush_iq_serv_bin_reg_3,
          obama_iq_serv_bin_reg_1,
          obama_iq_serv_bin_reg_2,
          obama_iq_serv_bin_reg_3,
          se = list(bush_iq_cas_bin_se_1,
                    bush_iq_cas_bin_se_2,
                    bush_iq_cas_bin_se_3,
                    bush_iq_serv_bin_se_1,
                    bush_iq_serv_bin_se_2,
                    bush_iq_serv_bin_se_3,
                    obama_iq_serv_bin_se_1,
                    obama_iq_serv_bin_se_2,
                    obama_iq_serv_bin_se_3),
          column.labels = rep(c("Combined", "Black", "White"), 3),
          covariate.labels = c("Black", "Casualty Sens.", "Military Family",
                               rep(NA, 23), "Constant"),
          omit = c("poll"),
          omit.labels = c("Poll FE"),
          add.lines = list(c("Polls", 18, 18, 18, 2,2,2,3,3,3)),
          omit.stat = c("adj.rsq","f", "ser"),
          out = "table_A16_cas_bin_regression_new.tex")

###################################################################
#
# 14) Costs of War - Spending (Bivariate)
#
###################################################################

###################################################################
# 14a) Regressions for Table A17, Models 1-3 (Spending, Bush)
###################################################################

# Model 1 - correlates of support for military spending, Bush
bush_iq_spend_bin_reg_1<-lm(dom_over_def ~ black+
                              factor(poll), data = cces_bush)

coeftest(bush_iq_spend_bin_reg_1, vcovHC(bush_iq_spend_bin_reg_1, 
                                         type = "HC1"))

bush_iq_spend_bin_se_1<-sqrt(diag(vcovHC(bush_iq_spend_bin_reg_1, 
                                         type = "HC1")))

# Model 2 - correlates of support for war, black sample
bush_iq_spend_bin_reg_2<-lm(iraq_support ~ dom_over_def+
                              factor(poll), data = cces_bush[which(cces_bush$black_aa==1),])

coeftest(bush_iq_spend_bin_reg_2, vcovHC(bush_iq_spend_bin_reg_2, 
                                         type = "HC1"))

bush_iq_spend_bin_se_2<-sqrt(diag(vcovHC(bush_iq_spend_bin_reg_2, 
                                         type = "HC1")))

# Model 3 - correlates of support for war, white sample
bush_iq_spend_bin_reg_3<-lm(iraq_support ~ dom_over_def+
                              factor(poll), data = cces_bush[which(cces_bush$black_aa==0),])

coeftest(bush_iq_spend_bin_reg_3, vcovHC(bush_iq_spend_bin_reg_3, 
                                         type = "HC1"))

bush_iq_spend_bin_se_3<-sqrt(diag(vcovHC(bush_iq_spend_bin_reg_3, 
                                         type = "HC1")))


# Mediation Analysis - Spending, Bush
set.seed(02138)

med.fit.bush.spend.bin <- lm(dom_over_def ~ black + factor(poll), 
                             data = cces_bush[-which(is.na(cces_bush$iraq_support)),])
out.fit.bush.spend.bin <- lm(iraq_support ~ black + dom_over_def + factor(poll),
                             data = cces_bush[-which(is.na(cces_bush$iraq_support)),])
med.out.bush.spend.bin <- mediate(med.fit.bush.spend.bin, 
                                  out.fit.bush.spend.bin, 
                                  treat = "black", mediator = "dom_over_def",
                                  robustSE = TRUE, sims = 250)
summary(med.out.bush.spend.bin)
plot(med.out.bush.spend.bin, xlim = c(-0.4, 0.1))

###################################################################
# 14b) Regressions for Table A17, Models 4-6 (Spending, Obama)
###################################################################

# Model 4 - correlates of support for military spending, Obama
obama_iq_spend_bin_reg_1<-lm(dom_over_def ~ black+
                               factor(poll), data = cces_obama)

coeftest(obama_iq_spend_bin_reg_1, vcovHC(obama_iq_spend_bin_reg_1, 
                                          type = "HC1"))

obama_iq_spend_bin_se_1<-sqrt(diag(vcovHC(obama_iq_spend_bin_reg_1, 
                                          type = "HC1")))

# Model 5 - correlates of support for war, black sample
obama_iq_spend_bin_reg_2<-lm(iraq_support ~ dom_over_def+
                               factor(poll), data = cces_obama[which(cces_obama$black_aa==1),])

coeftest(obama_iq_spend_bin_reg_2, vcovHC(obama_iq_spend_bin_reg_2, 
                                          type = "HC1"))

obama_iq_spend_bin_se_2<-sqrt(diag(vcovHC(obama_iq_spend_bin_reg_2, 
                                          type = "HC1")))

# Model 6 - correlates of support for war, white sample
obama_iq_spend_bin_reg_3<-lm(iraq_support ~ dom_over_def+
                               factor(poll), 
                             data = cces_obama[which(cces_obama$black_aa==0),])

coeftest(obama_iq_spend_bin_reg_3, vcovHC(obama_iq_spend_bin_reg_3, 
                                          type = "HC1"))

obama_iq_spend_bin_se_3<-sqrt(diag(vcovHC(obama_iq_spend_bin_reg_3, 
                                          type = "HC1")))

# values
coef(bush_iq_spend_bin_reg_1)[2]+1.96*bush_iq_spend_bin_se_1[2]
coef(bush_iq_spend_bin_reg_1)[2]-1.96*bush_iq_spend_bin_se_1[2]

coef(obama_iq_spend_bin_reg_1)[2]+1.96*obama_iq_spend_bin_se_1[2]
coef(obama_iq_spend_bin_reg_1)[2]-1.96*obama_iq_spend_bin_se_1[2]

# Mediation Analysis - Spending, Obama
set.seed(02138)

med.fit.obama.spend.bin <- lm(dom_over_def ~ black + factor(poll), 
                              data = cces_obama[-which(is.na(cces_obama$iraq_support)),])
out.fit.obama.spend.bin <- lm(iraq_support ~ black + dom_over_def + factor(poll), data = cces_obama[-which(is.na(cces_obama$iraq_support)),])
med.out.obama.spend.bin <- mediate(med.fit.obama.spend.bin, out.fit.obama.spend.bin, treat = "black", mediator = "dom_over_def",
                                   robustSE = TRUE, sims = 250)
summary(med.out.obama.spend.bin)
plot(med.out.obama.spend.bin, xlim = c(-0.4, 0.1))

###################################################################
# 14c) Generates Table A17, Models 1-6 (Spending & War Bin)
###################################################################
stargazer(bush_iq_spend_bin_reg_1,
          bush_iq_spend_bin_reg_2,
          bush_iq_spend_bin_reg_3,
          obama_iq_spend_bin_reg_1,
          obama_iq_spend_bin_reg_2,
          obama_iq_spend_bin_reg_3,
          se = list(bush_iq_spend_bin_se_1,
                    bush_iq_spend_bin_se_2,
                    bush_iq_spend_bin_se_3,
                    obama_iq_spend_bin_se_1,
                    obama_iq_spend_bin_se_2,
                    obama_iq_spend_bin_se_3),
          covariate.labels = c("Black", "Dom. Spending",
                               rep(NA, 3), "Constant"),
          column.labels = rep(c("Combined", "Black", "White"),2),
          omit = c("poll"),
          omit.labels = c("Poll FE"),
          add.lines = list(c("Polls", 2,2,2,3,3,3)),
          omit.stat = c("adj.rsq","f", "ser"),
          out = "table_A17_spending_bin_regression_new.tex")

###################################################################
# 13d) Figure A2, Mediation Plot for Spending & War
###################################################################
pdf("spending_bin_mediation_plots_new.pdf", width = 4, height = 4)
par(mar=c(2, 4, 4, 2),
    mfcol = c(2, 1))
plot(med.out.bush.spend.bin, xlim = c(-0.4, 0.1), cex.axis = .8)
title(main = list("Spending (Bush)", cex = 1,
                  col = "black", font = 1), line = 1, adj = 0)
plot(med.out.obama.spend.bin, xlim = c(-0.4, 0.1), cex.axis = .8)
title(main = list("Spending (Obama)", cex = 1,
                  col = "black", font = 1), line = 1, adj = 0)

dev.off()

###################################################################
#
# 15) Prerequisites for War - Alienation (Binary)
#
###################################################################

#############################################################################
# 15a) Regressions for Figure A3 & Table A18, Models 1-3 (Alienation, Bush)
#############################################################################

# Model 1 - correlates of alienation, Bush
Bush1alien_IQ_bin <- lm(alienation_binary ~ 
                          black+factor(poll), data = Bush_data_IQ_BW)
Bush1alien_IQ_bin_se<-sqrt(diag(vcovHC(Bush1alien_IQ_bin, type = "HC1")))

# Model 2 - correlates of support for war, black sample
Bush2alien_IQ_bin <- lm(outcome ~ alienation_binary+
                          factor(poll), data = Bush_data_IQ_BW[which(Bush_data_IQ_BW$black==1),])
Bush2alien_IQ_bin_se<-sqrt(diag(vcovHC(Bush2alien_IQ_bin, type = "HC1")))

# Model 3 - correlates of support for war, white sample
Bush3alien_IQ_bin <- lm(outcome ~ alienation_binary+
                          factor(poll), data = Bush_data_IQ_BW[which(Bush_data_IQ_BW$black==0),])
Bush3alien_IQ_bin_se<-sqrt(diag(vcovHC(Bush3alien_IQ_bin, type = "HC1")))

coef(Bush1alien_IQ_bin)[2]+1.96*Bush1alien_IQ_bin_se[2]
coef(Bush1alien_IQ_bin)[2]-1.96*Bush1alien_IQ_bin_se[2]

# Mediation - Bush & Alienation
set.seed(02138)

med.fit.bush.alien.bin <- lm(alienation_binary ~ black + 
                               factor(poll), data = Bush_data_IQ_BW[!is.na(Bush_data_IQ_BW$outcome),])
out.fit.bush.alien.bin <- lm(outcome ~ black + alienation_binary + 
                               factor(poll), data = Bush_data_IQ_BW[!is.na(Bush_data_IQ_BW$outcome),])
med.out.bush.alien.bin <- mediate(med.fit.bush.alien.bin, out.fit.bush.alien.bin, 
                                  treat = "black", mediator = "alienation_binary",
                                  robustSE = TRUE, sims = 250)
summary(med.out.bush.alien.bin)
plot(med.out.bush.alien.bin, xlim = c(-0.4, 0.1))


##############################################################################
# 15b) Regressions for Figure 3A3 & Table A18, Models 4-6 (Alienation, Obama)
##############################################################################

# Model 4 - correlates of alienation, Obama
Obama1alien_IQ_bin <- lm(alienation_binary ~ 
                           black+factor(poll), data = Obama_data_IQ_BW)
Obama1alien_IQ_bin_se<-sqrt(diag(vcovHC(Obama1alien_IQ_bin, type = "HC1")))

# Model 5 - correlates of support for war, black sample
Obama2alien_IQ_bin <- lm(outcome ~ alienation_binary +
                           +factor(poll), data = Obama_data_IQ_BW[which(Obama_data_IQ_BW$black==1),])
Obama2alien_IQ_bin_se<-sqrt(diag(vcovHC(Obama2alien_IQ_bin, type = "HC1")))

# Model 6 - correlates of support for war, white sample
Obama3alien_IQ_bin <- lm(outcome ~ alienation_binary +
                           factor(poll), data = Obama_data_IQ_BW[which(Obama_data_IQ_BW$black==0),])
Obama3alien_IQ_bin_se<-sqrt(diag(vcovHC(Obama3alien_IQ_bin, type = "HC1")))

summary(Obama3alien_IQ_bin)

#
coef(Obama1alien_IQ_bin)[2]+1.96*Obama1alien_IQ_bin_se[2]
coef(Obama1alien_IQ_bin)[2]-1.96*Obama1alien_IQ_bin_se[2]

# Mediation - Obama & Alienation
set.seed(02138)

med.fit.obama.alien.bin <- lm(alienation_binary ~ black + 
                                factor(poll), 
                              data = Obama_data_IQ_BW[which(!is.na(Obama_data_IQ_BW$outcome)),])
out.fit.obama.alien.bin <- lm(outcome ~ black + alienation_binary + 
                                factor(poll), 
                              data = Obama_data_IQ_BW[which(!is.na(Obama_data_IQ_BW$outcome)),])
med.out.obama.alien.bin <- mediate(med.fit.obama.alien.bin, out.fit.obama.alien.bin, 
                                   treat = "black", mediator = "alienation_binary",
                                   robustSE = TRUE, sims = 250)
summary(med.out.obama.alien.bin)
plot(med.out.obama.alien.bin, xlim = c(-0.4, 0.1))

###################################################################
# 15c) Generates Table A18, Mediation Plot for Spending & War Bin
###################################################################
stargazer(Bush1alien_IQ_bin,Bush2alien_IQ_bin,Bush3alien_IQ_bin,
          Obama1alien_IQ_bin,Obama2alien_IQ_bin,Obama3alien_IQ_bin,
          se = list(Bush1alien_IQ_bin_se,Bush2alien_IQ_bin_se,Bush3alien_IQ_bin_se,
                    Obama1alien_IQ_bin_se,Obama2alien_IQ_bin_se,Obama3alien_IQ_bin_se),
          covariate.labels = c("Black", "Alienation",
                               rep(NA, 43), "Constant"),
          column.labels = rep(c("Combined", "Black", "White"),2),
          omit = c("poll"),
          omit.labels = c("Poll FE"),
          add.lines = list(c("Polls", 39,39,39,4,4,4)),
          omit.stat = c("adj.rsq","f", "ser"),
          out = "table_A18_alienation_bin_regression_new.tex")

###################################################################
# 15d) Figure A3, Mediation Plot for Alienation & War
###################################################################
pdf("alien_bin_mediation_plots_new.pdf", width = 4, height = 4)
par(mar=c(2, 4, 4, 2),
    mfcol = c(2, 1))
plot(med.out.bush.alien.bin, xlim = c(-0.4, 0.1), cex.axis = .8)
title(main = list("Alienation (Bush)", cex = 1,
                  col = "black", font = 1), line = 1, adj = 0)
plot(med.out.obama.alien.bin, xlim = c(-0.4, 0.1), cex.axis = .8)
title(main = list("Alienation (Obama)", cex = 1,
                  col = "black", font = 1), line = 1, adj = 0)
dev.off()

########################################################
#
# 18) Sensitivity Analysis
#
########################################################

# Sensitivity - Bush Casualties Iraq

set.seed(02138)
sens.out.bush.cas <- medsens(med.out.bush.cas, rho.by = 0.1, effect.type = "indirect", sims = 100)
summary(sens.out.bush.cas)
plot(sens.out.bush.cas)

set.seed(02138)
sens.out.bush.cas.bin <- medsens(med.out.bush.cas.bin, rho.by = 0.1, effect.type = "indirect", sims = 100)
summary(sens.out.bush.cas.bin)
plot(sens.out.bush.cas.bin)

# Sensitivity - Bush Serving Iraq

set.seed(02138)
sens.out.bush.serv <- medsens(med.out.bush.serv, rho.by = 0.1, effect.type = "indirect", sims = 100)
summary(sens.out.bush.serv)
plot(sens.out.bush.serv)

set.seed(02138)
sens.out.bush.serv.bin <- medsens(med.out.bush.serv.bin, rho.by = 0.1, effect.type = "indirect", sims = 100)
summary(sens.out.bush.serv.bin)
plot(sens.out.bush.serv.bin)

# Sensitivity - Obama Serving Iraq

set.seed(02138)
sens.out.obama.serv <- medsens(med.out.obama.serv, rho.by = 0.1, effect.type = "indirect", sims = 100)
summary(sens.out.obama.serv)
plot(sens.out.obama.serv)

set.seed(02138)
sens.out.obama.serv.bin <- medsens(med.out.obama.serv.bin, rho.by = 0.1, effect.type = "indirect", sims = 100)
summary(sens.out.obama.serv.bin)
plot(sens.out.obama.serv.bin)


# Sensitivity - Bush Spending Iraq

set.seed(02138)
sens.out.bush.spend <- medsens(med.out.bush.spend, rho.by = 0.1, effect.type = "indirect", sims = 100)
summary(sens.out.bush.spend)
plot(sens.out.bush.spend)

set.seed(02138)
sens.out.bush.spend.bin <- medsens(med.out.bush.spend.bin, rho.by = 0.1, effect.type = "indirect", sims = 100)
summary(sens.out.bush.spend.bin)
plot(sens.out.bush.spend.bin)

# Sensitivity - Obama Spending Iraq

set.seed(02138)
sens.out.obama.spend <- medsens(med.out.obama.spend, rho.by = 0.1, effect.type = "indirect", sims = 100)
summary(sens.out.obama.spend)
plot(sens.out.obama.spend)

set.seed(02138)
sens.out.obama.spend.bin <- medsens(med.out.obama.spend.bin, rho.by = 0.1, effect.type = "indirect", sims = 100)
summary(sens.out.obama.spend.bin)
plot(sens.out.obama.spend.bin)


# Sensitivity - Bush Alienation Iraq

set.seed(02138)
sens.out.bush.alien <- medsens(med.out.bush.alien, rho.by = 0.1, effect.type = "indirect", sims = 100)
summary(sens.out.bush.alien)
plot(sens.out.bush.alien)

set.seed(02138)
sens.out.bush.alien.bin <- medsens(med.out.bush.alien.bin, rho.by = 0.1, effect.type = "indirect", sims = 100)
summary(sens.out.bush.alien.bin)
plot(sens.out.bush.alien.bin)

# Sensitivity - Obama Alienation Iraq

set.seed(02138)
sens.out.obama.alien <- medsens(med.out.obama.alien, rho.by = 0.1, effect.type = "indirect", sims = 100)
summary(sens.out.obama.alien)
plot(sens.out.obama.alien)

set.seed(02138)
sens.out.obama.alien.bin <- medsens(med.out.obama.alien.bin, rho.by = 0.1, effect.type = "indirect", sims = 100)
summary(sens.out.obama.alien.bin)
plot(sens.out.obama.alien.bin)

########################################################
#
# 17) Specific statistics
#
########################################################

# pg. 15 --> n of ABC polls
# casualties - Iraq - Bush
nrow(abc_iraq_combine)

table(abc_iraq_combine$black)
table(abc_iraq_combine$black)/nrow(abc_iraq_combine)

# consistency of casualty support 
casualty_consist<-table(abc_iraq_combine$casual_sens, abc_iraq_combine$support_war)/sum(table(abc_iraq_combine$casual_sens, abc_iraq_combine$support_war))
casualty_consist
casualty_consist[2,1]+casualty_consist[1,2]

# avg. size of polls
mean(table(abc_iraq_combine$poll_code))

# min and max
min(table(abc_iraq_combine$poll_code))
max(table(abc_iraq_combine$poll_code))

# casualties - Afghanistan - Bush
nrow(poll_afgh)
table(poll_afgh$black)
table(poll_afgh$black)/nrow(poll_afgh)

# descriptive stats for alienation polls
# alienation - Iraq - Bush

nrow(Bush_data_IQ_BW)
mean(Bush_data_IQ_BW$black)

# alienation - Iraq - Obama

nrow(Obama_data_IQ_BW)
mean(Obama_data_IQ_BW$black)

# descriptive stats for cces polls
# CCES - Bush

nrow(cces_bush)
mean(cces_bush$black)

# CCES - Obama

nrow(cces_obama)
mean(cces_obama$black)

########################################################
#
# 18) Gulf War barplot
#
########################################################

gulf_burden<-read.dta("black_burden_gulf_war.dta")

table(gulf_burden$Q861J)

gulf_burden <- gulf_burden[gulf_burden$Q861J == "white" | gulf_burden$Q861J == "black",]

table(gulf_burden$Q861J)

gulf_burden$black_burden <- ifelse(gulf_burden$Q9D == "MODERATELY AGREE" | gulf_burden$Q9D == "STRONGLY AGREE", 1, 0)
gulf_burden$black <- ifelse(gulf_burden$Q861J == "black", 1, 0)

table(gulf_burden$black_burden, gulf_burden$black)
table(gulf_burden$black)

b_burden_agree <- 120/206 
w_burden_agree <- 175/741


barplot(c(b_burden_agree, w_burden_agree), col = c("black", "blue"))

# Create data
burden_data <- data.frame(
  name=c("White", "Black/Af.-Am.") ,  
  value=c(w_burden_agree, b_burden_agree)
)

# Barplot
pdf("burden_garamond.pdf", width=6, height=4)
burden_plot <- ggplot(burden_data, aes(x=name, y=value, fill = as.factor(name))) + 
  ggtitle("Agree: Black Americans will bear an unfair burden in the war.")+
  geom_bar(stat="identity")+theme_bw()+
  xlab("") + 
  ylab("% of Respondents") + 
  scale_fill_manual(values = c("black", "blue"))+
  scale_y_continuous("% of Respondents", limits = c(0,.6),
                     breaks = seq(0, .6, by = 0.2),
                     labels = c("0%", "20%", "40%", "60%"))

burden_plot +theme_set(theme_bw(base_size = 10, base_family = "EB Garamond" )) + theme(legend.position = "none")
dev.off()

setEPS()                                             
postscript("burden_garamond.eps", width=6, height=4)
burden_plot <- ggplot(burden_data, aes(x=name, y=value, fill = as.factor(name))) + 
  ggtitle("Agree: Black Americans will bear an unfair burden in the war.")+
  geom_bar(stat="identity")+theme_bw()+
  xlab("") + 
  ylab("% of Respondents") + 
  scale_fill_manual(values = c("black", "blue"))+
  scale_y_continuous("% of Respondents", limits = c(0,.6),
                     breaks = seq(0, .6, by = 0.2),
                     labels = c("0%", "20%", "40%", "60%"))

burden_plot +theme_set(theme_bw(base_size = 10, base_family = "EB Garamond" )) + theme(legend.position = "none")
dev.off()

